knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE, fig.align='center')

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.4.4
library(plyr)
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.4.4
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(plotrix)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:plyr':
## 
##     here
## The following object is masked from 'package:base':
## 
##     date
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(RCurl)
## Warning: package 'RCurl' was built under R version 3.4.4
## Loading required package: bitops
library(chron)
## 
## Attaching package: 'chron'
## The following objects are masked from 'package:lubridate':
## 
##     days, hours, minutes, seconds, years
library(lubridate)
library(lme4)
## Loading required package: Matrix
library(emmeans)
## Warning: package 'emmeans' was built under R version 3.4.4
## NOTE: As of emmeans versions > 1.2.3,
##       The 'cld' function will be deprecated in favor of 'CLD'.
##       You may use 'cld' only if you have package:multcomp attached.
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## 
## Attaching package: 'multcomp'
## The following object is masked from 'package:emmeans':
## 
##     cld
library(devtools)
## Warning: package 'devtools' was built under R version 3.4.4
## 
## Attaching package: 'devtools'
## The following object is masked from 'package:emmeans':
## 
##     test
library(logihist)
library(popbio)
## Warning: package 'popbio' was built under R version 3.4.4
library(car)
## 
## Attaching package: 'car'
## The following objects are masked from 'package:carData':
## 
##     Guyer, UN, Vocab
## The following object is masked from 'package:dplyr':
## 
##     recode
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:emmeans':
## 
##     lsmeans
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
## 
## Attaching package: 'sjPlot'
## The following objects are masked from 'package:cowplot':
## 
##     plot_grid, save_plot
library(sjmisc)

Project background

Environmental

Sea level correction

Corals were collected across in summer and winter 2016 across several different days with each period. To correct for differences in the depth of coral colonies we used NOAA tide data to correct all depth measurements to Mean Sea Level (MSL). This approach was done in Innis et al. 2018 on a larger collection of the coral samples–which have been subset in the current study.

#########################
# Sea level correction
#########################

#### attach data files
data<-read.csv("data/mastersheet_PanKBAY.csv") # master file
qPCR.Innis<-read.csv("data/PanKBay_summer_qPCR.csv") # qPCR from summer (Innis et al. 2018)

JuneTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160801-20160831.csv")
DecemberTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20161201-20161222.csv")

Tide<-rbind(JuneTide, JulyTide, AugustTide, DecemberTide) # all MSL in meters

Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
# use attr(as.POSIXlt(Tide$Time),"tzone") to confirm in PST

Tide<-Tide[, c(-5:-8)] #remove unnecessary columns

data$date.time<- as.POSIXct(paste(as.character(data$Date), as.character(data$Time.of.collection)),
                 format="%m/%d/%y %H:%M", tz="Pacific/Honolulu") # make sure time in HST for master

data$Time.r<-round_date(data$date.time,unit="6 minutes")
Tide$Time.r <- Tide$Time
data<-merge(data, Tide, by="Time.r", all.x=T)
data$newDepth <- data$Depth..m-data$TideHT # new depth in METERS

Light parameters

DLI 2m all sites

This figure shows the daily integrated light (DLI mol photon m-2 d-1) at long-term light loggers at 2m depth near the locations where corals were collected in northern and southern sectors of Kāne’ohe Bay near shore (east) further off shore (west).

# DLI
##### all DLI for graphs
files <- list.files(path="data/environmental/temp and light/Jun_DecPAR/all DLI", pattern = "csv$", full.names = T)
tables <- lapply(files, read.csv, header = TRUE)
DLI<-do.call(rbind, tables)
DLI=DLI[-1] # remove junk column
DLI$Date<-as.Date(DLI$Date) # fix date
DLI<-DLI[order(DLI$Date),] # order by Date

# make a date sequence for entire study
all.date<-as.data.frame(seq(as.Date("2016-06-10"), as.Date("2017-01-12"), "days"))
colnames(all.date)="Date"

# merge the DLI data and the date sequence to make a complete df through time
DLI.3site<-merge(all.date, DLI, by="Date", all.x=T)
R10<-read.csv("data/environmental/temp and light/Jun_DecPAR/all DLI/Reef 10/DLI.Rf102016.csv")
R10<-R10[-1]; R10$Date<-as.Date(R10$Date)

DLI.4site<-merge(DLI.3site, R10, by="Date", all.x=T)
DLI.4site$month <- months(as.Date(DLI.4site$Date)) # makes a month column
DLI.4site<-DLI.4site[, c(1,6, 2:5)]

# write.csv(DLI.4site, "data/environmental/temp and light/Jun_DecPAR/All.DLI.csv")

##### Figure
All.DLI<-DLI.4site
reefcols=c("mediumseagreen", "dodgerblue", "salmon", "orchid")
par(mar=c(2,3.6,1,1.3), mgp=c(2,0.5,0))
layout(matrix(c(1,2,3,4), 2,2, byrow=TRUE))

plot(Rf44.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[1])
legend("topright", lty=1, col=reefcols[1], legend="Northwest", lwd=1.5, bty="n", cex=0.8, 
       x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")

plot(Rf42.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[2])
legend("topright", lty=1, col=reefcols[2], legend="Northeast", lwd=1.5, bty="n", cex=0.8, 
       x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")

plot(Rf10.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[3])
legend("topright", lty=1, col=reefcols[3], legend="Southwest", lwd=1.5, bty="n", cex=0.8, 
       x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")

plot(HIMB.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[4])
legend("topright", lty=1, col=reefcols[4], legend="HIMB", lwd=1.5, bty="n", cex=0.8, 
       x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
**Figure** *DLI at the four reef areas from summer to winter 2016 at 2m depths. Data at 2m was recorded from loggers*

Figure DLI at the four reef areas from summer to winter 2016 at 2m depths. Data at 2m was recorded from loggers

dev.copy(pdf, "figures/environmental/All.DLI.pdf", width=6.5, height=4)
dev.off()

Model results testing differences in DLI at 2 m depth at all sites.

mod.light<-DLI.4site
Rf44DLI.mod<-mod.light[1:3]; Rf44DLI.mod$Location<-"Rf44"; 
    colnames(Rf44DLI.mod)<-c("Date", "month", "DLI", "Location")
Rf42DLI.mod<-mod.light[,c(1:2,4)]; Rf42DLI.mod$Location<-"Rf42"; 
    colnames(Rf42DLI.mod)<-c("Date", "month", "DLI", "Location")
HIMBDLI.mod<-mod.light[,c(1:2,5)]; HIMBDLI.mod$Location<-"HIMB"; 
    colnames(HIMBDLI.mod)<-c("Date", "month", "DLI", "Location")
Rf10DLI.mod<-mod.light[,c(1:2,6)]; Rf10DLI.mod$Location<-"Rf10"; 
    colnames(Rf10DLI.mod)<-c("Date", "month", "DLI", "Location")

All.DLI.mod<-rbind(Rf44DLI.mod, Rf42DLI.mod, HIMBDLI.mod, Rf10DLI.mod)
All.DLI.mod$Location<-as.factor(All.DLI.mod$Location)
All.DLI.mod$month<-as.factor(All.DLI.mod$month)
All.DLI.mod$Season<-ifelse(All.DLI.mod$month== "June" | 
                             All.DLI.mod$month== "July" |
                                All.DLI.mod$month== "August" |
                                  All.DLI.mod$month== "September", "dry season", "wet season")
All.DLI.mod$Season<-as.factor(All.DLI.mod$Season)

mod<-lmer(DLI~Location*Season + (1|Date), data=All.DLI.mod); anova(mod, type=2)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##                 Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Location        4378.8 1459.58     3 529.58 134.674 < 2.2e-16 ***
## Season          1040.9 1040.91     1 210.48  96.043 < 2.2e-16 ***
## Location:Season  490.9  163.62     3 530.91  15.097 1.924e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(mod))
posthoc<-emmeans(mod, ~Location|Season)
CLD(posthoc, Letters=letters)
## Season = dry season:
##  Location    emmean        SE     df  lower.CL  upper.CL .group
##  Rf10      8.149681 0.4014748 504.55  7.360912  8.938449  a    
##  Rf44     10.829854 0.4450350 590.91  9.955811 11.703896   b   
##  HIMB     14.748648 0.4016675 503.81 13.959498 15.537797    c  
##  Rf42     14.997825 0.4601858 615.62 14.094100 15.901549    c  
## 
## Season = wet season:
##  Location    emmean        SE     df  lower.CL  upper.CL .group
##  Rf10      3.901772 0.4473755 539.33  3.022959  4.780584  a    
##  HIMB      8.607241 0.4662620 551.22  7.691373  9.523109   b   
##  Rf44      9.128302 0.4662620 551.22  8.212435 10.044170   b   
##  Rf42      9.448888 0.4707121 559.55  8.524310 10.373467   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-emmeans(mod, ~Season)
CLD(posthoc, Letters=letters)
##  Season        emmean        SE     df  lower.CL  upper.CL .group
##  wet season  7.771551 0.3348801 206.76  7.111333  8.431768  a    
##  dry season 12.181502 0.3071824 206.48 11.575885 12.787118   b   
## 
## Results are averaged over the levels of: Location 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
**Figure** DLI model output

Figure DLI model output

PAR and DLI at 3 depths

With Depth as a factor (<1m, 2m, 8m), we can produce a continuous graph of light (PAR or DLI) at 2m and use model intercepts to predict the light at the shallow and deep zones (ca. <1m and 8m, respectfully). Here, the model is \(lm(log.PAR\) ~ \(Depth.factor)\).

# Plot it!!
# Reef 44
par(mfrow=c(2,2), mar=c(4,5,2,2))

plot(y=All.PAR.df$Rf44.shall, x=All.PAR.df$timestamp, type="l", col="palegreen3", 
     cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
     , xlab="", main="Northwest", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf44.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf44.deep, x=All.PAR.df$timestamp, type="l", col="palegreen4"))
legend("topright", lty=1, col=c("palegreen3", "gray", "palegreen4"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)


# Plot it!!
# Reef 42
plot(y=All.PAR.df$Rf42.shall, x=All.PAR.df$timestamp, type="l", col="mediumturquoise", 
     cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
     , xlab="", main="Northeast", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf42.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf42.deep, x=All.PAR.df$timestamp, type="l", col="dodgerblue2"))
legend("topright", lty=1, col=c("mediumturquoise", "gray", "dodgerblue2"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)

# Plot it!!
# Reef 10
plot(y=All.PAR.df$Rf10.shall, x=All.PAR.df$timestamp, type="l", col="orangered", 
     cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
     , xlab="Dates", main="Southwest", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf10.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf10.deep, x=All.PAR.df$timestamp, type="l", col="lightsalmon3"))
legend("topright", lty=1, col=c("orangered", "gray", "lightsalmon3"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)

# Plot it!
# HIMB
plot(y=All.PAR.df$HIMB.shall, x=All.PAR.df$timestamp, type="l", col="mediumorchid2", 
     cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
     , xlab="Dates", main="HIMB", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$HIMB.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$HIMB.deep, x=All.PAR.df$timestamp, type="l", col="plum4"))
legend("topright", lty=1, col=c("mediumorchid2", "gray", "plum4"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
**Figure** *PAR at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations*

Figure PAR at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations

dev.copy(pdf, "figures/environmental/PAR.all depths.pdf", height=7, width=8)
dev.off()
#######
####### calculate DLI and compare
#######
df<-All.PAR.df
df$timestamp<-strptime(df$timestamp, format="%Y-%m-%d %H:%M:%S")
df$timestamp<-as.Date(df$timestamp, format="%Y-%m-%d") # change to DATE ONLY format to calulate range, means
df[is.na(df)] <- 0


df.split <- split(df, f=df$timestamp
                  < as.Date("2016-06-10", format="%Y-%m-%d")) # split df by date

df.dli<-aggregate(data.frame(Rf42.shall.DLI=df.split[[1]]$Rf42.shall*0.0864,
                             Rf42.mid.DLI=df.split[[1]]$Rf42.mid*0.0864,
                             Rf42.deep.DLI=df.split[[1]]$Rf42.deep*0.0864,
                             Rf44.shall.DLI=df.split[[1]]$Rf44.shall*0.0864,
                             Rf44.mid.DLI=df.split[[1]]$Rf44.mid*0.0864,
                             Rf44.deep.DLI=df.split[[1]]$Rf44.deep*0.0864,
                             HIMB.shall.DLI=df.split[[1]]$HIMB.shall*0.0864,
                             HIMB.mid.DLI=df.split[[1]]$HIMB.mid*0.0864,
                             HIMB.deep.DLI=df.split[[1]]$HIMB.deep*0.0864,
                             Rf10.shall.DLI=df.split[[1]]$Rf10.shall*0.0864,
                             Rf10.mid.DLI=df.split[[1]]$Rf10.mid*0.0864,
                             Rf10.deep.DLI=df.split[[1]]$Rf10.deep*0.0864),
                  by=list(Date=df.split[[1]]$timestamp), FUN=mean)


# make a date sequence for entire study
all.date.time<-as.data.frame(seq(
  from=as.POSIXct("2016-06-10", tz="HST"),
  to=as.POSIXct("2017-01-12", tz="HST"),
  by="1 d"))  
colnames(all.date.time)[1]<-"Date"
all.date.time$Date<-strptime(all.date.time$Date, format="%Y-%m-%d")
all.date.time$Date<-as.Date(all.date.time$Date, format="%Y-%m-%d")

df.dli<-merge(all.date.time, df.dli, by="Date", all.x=T)
df.dli[df.dli<=0] <- NA

#####################
# Plot it!!
# Reef 44
par(mfrow=c(2,2), mar=c(4,5,2,2))

plot(y=df.dli$Rf44.shall, x=df.dli$Date, type="h", col="palegreen3", 
     cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
     , xlab="", main="Northwest", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf44.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf44.deep, x=df.dli$Date, type="h", col="palegreen4"))
legend("topright", lty=1, col=c("palegreen3", "gray", "palegreen4"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)


# Plot it!!
# Reef 42
plot(y=df.dli$Rf42.shall, x=df.dli$Date, type="h", col="mediumturquoise", 
     cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
     , xlab="", main="Northeast", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf42.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf42.deep, x=df.dli$Date, type="h", col="dodgerblue2"))
legend("topright", lty=1, col=c("mediumturquoise", "gray", "dodgerblue2"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)

# Plot it!!
# Reef 10
plot(y=df.dli$Rf10.shall, x=df.dli$Date, type="h", col="orangered", 
     cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
     , xlab="Dates", main="Southwest", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf10.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf10.deep, x=df.dli$Date, type="h", col="lightsalmon3"))
legend("topright", lty=1, col=c("orangered", "gray", "lightsalmon3"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)

# Plot it!
# HIMB
plot(y=df.dli$HIMB.shall, x=df.dli$Date, type="h", col="mediumorchid2", 
     cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
     , xlab="Dates", main="Southeast", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$HIMB.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$HIMB.deep, x=df.dli$Date, type="h", col="plum4"))
legend("topright", lty=1, col=c("mediumorchid2", "gray", "plum4"), legend=c("<1m", "2m", "8m"), 
       lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
**Figure** *DLI at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations*

Figure DLI at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations

dev.copy(pdf, "figures/environmental/DLIcalc.all depths.pdf", height=7, width=8)
dev.off()

Attenuation and Light-at-Depth

A periodic deployment of loggers (October and November) at ca. <1m and ca. 7m were used to model the relationship between 2m light and light at other depths.

Using the logger data, the difference in light at 1 and 7m was calculated relative to 2m (i.e, PAR baseline (2m) - PAR other), then the difference in depth at 1 and 7m was calculated realtive to logger depth i.e, depth logger other - depth logger (2m baseline). A no-intercept model was used to fit difference in log(DLI) and depth (i.e, log(delta-DLI) and delta-depth) with depth as a continous/numeric variable to determine kd: \(lm(delta.log.DLI\)~\(delta.Depth + 0)\).

Site-specific kd values were used to generate a continuous variable of “light at depth”. At each time period (summer or winter) the daily light integral was averaged over 3 months: summer 2016 (June, July, August) and winter 2016 (November, December, January). Discrete monthly seasonal-mean DLI means was used to generate an estimate of light at depth for corals at each site, as this integrates the site and seasoanl variance.

Light at Depth was calculated following Beer-Lambert: \(Ez{_d} = Ez_{_0}e^{(-k_d *d)}\), where \(d\) is depth, \(k_d\) is the attenuation coefficient for each site, \(Ez{_d}\) is light at depth \(d\) and \(Ez{_0}\) is measured light. Calculations were relative to light at 2m depth (i.e., \(Ez_{_0}\) was the light at 2m depth) and \(d\) was the difference between coral depth relative to depth of the logger (i.e., 2m). light reaching depth d and Ez(0) is the incident light at the surface.

########################## # LIGHT AT DEPTH
#############
###### 
# using new depth and the site-specific kd (light attenuation) determine approximate "light at depth" for each colony sample. 

######

### attach necessary files
logger.depths<-read.csv("data/environmental/temp and light/light.logger.depths.csv") # depths for loggers
kd.all<-read.csv("data/environmental/kd.all.csv") # kds for each site

# Seasonal DLI used for "period of collection" light levels
month.DLI<-read.csv("data/environmental/temp and light/Jun_DecPAR/All.DLI_long.csv")

# corals collected in June-July-August use summer time DLI for these months as indicator of average DLI
summer.DLI<-month.DLI[(month.DLI$Month=="June" | month.DLI$Month=="July" | month.DLI$Month=="August"),]
winter.DLI<-month.DLI[(month.DLI$Month=="November" | month.DLI$Month=="December" |month.DLI$Month=="January"),]

# summer mean and SE dataframe
sum.mean<-aggregate(DLI~Site, summer.DLI, mean)
sum.SE<-aggregate(DLI~Site, summer.DLI, std.error)
sum.light.df<-cbind(sum.mean, sum.SE[2])
colnames(sum.light.df)=c("Site", "mean.DLI", "SE")
sum.light.df$Season<-as.factor("summer")

# winter mean and SE dataframe
wint.mean<-aggregate(DLI~Site, winter.DLI, mean)
wint.SE<-aggregate(DLI~Site, winter.DLI, std.error)
wint.light.df<-cbind(wint.mean, wint.SE[2])
colnames(wint.light.df)=c("Site", "mean.DLI", "SE")
wint.light.df$Season<-as.factor("winter")

season.DLI<-rbind(sum.light.df[,c(4,1,2:3)], wint.light.df[,c(4,1,2:3)]) # compiled means for DLI at ~2m

write.csv(season.DLI, "data/environmental/season.DLI.csv")
#######################################
### make new dataframe for calculations
df.light<- data[, c("Season", "Location", "newDepth")]

# make a column of "depth differences" relative to where ~2m logger was deployed
df.light$depth.diff<-ifelse(df.light$Location=="F1-R46", df.light$newDepth-1.8288,
                          ifelse(df.light$Location=="F8-R10", df.light$newDepth-1.8288, 
                                 ifelse(df.light$Location=="HIMB", df.light$newDepth-1.8288, 
                                        df.light$newDepth-2.4384))) # last statement for remaining site (Rf42)

# make a column for sample-specific light at depth (estimate) based on kd
# follow: #with 2m as baseline#  E(depth) = E(2m)*exp(-k_d*(depth - 2m))
# so that:     light at depth = DLI at mid.depth * exp (-kd *(delta shallow-deep in m))


df.light$Light<-ifelse(df.light$Location=="HIMB" & df.light$Season=="summer", 
                          season.DLI$mean.DLI[1]*exp(-kd.all$kd[1]*df.light$depth.diff), # summer HIMB
                   ifelse(df.light$Location=="F8-R10" & df.light$Season=="summer", 
                          season.DLI$mean.DLI[2]*exp(-kd.all$kd[2]*df.light$depth.diff), # summer R10
                   ifelse(df.light$Location=="R42" & df.light$Season=="summer", 
                          season.DLI$mean.DLI[3]*exp(-kd.all$kd[3]*df.light$depth.diff), # summer R42
                   ifelse(df.light$Location=="F1-R46" & df.light$Season=="summer", 
                                 season.DLI$mean.DLI[4]*exp(-kd.all$kd[4]*df.light$depth.diff), # summer R46
                          
                   ifelse(df.light$Location=="HIMB" & df.light$Season=="winter", 
                          season.DLI$mean.DLI[5]*exp(-kd.all$kd[1]*df.light$depth.diff), # winter HIMB
                   ifelse(df.light$Location=="F8-R10" & df.light$Season=="winter", 
                          season.DLI$mean.DLI[6]*exp(-kd.all$kd[2]*df.light$depth.diff), # winter R10
                   ifelse(df.light$Location=="R42" & df.light$Season=="winter", 
                          season.DLI$mean.DLI[7]*exp(-kd.all$kd[3]*df.light$depth.diff), # winter R42
                          season.DLI$mean.DLI[8]*exp(-kd.all$kd[4]*df.light$depth.diff)) # winter R46
                            ))))))

###### plot of light x depth by season
df.light$Location <- factor(df.light$Location, levels = c("F1-R46", "R42", "F8-R10", "HIMB"))

plot.by.sites=ggplot(df.light, aes(Light)) + geom_density(aes(fill=Location), alpha=0.3, position = 'stack') + scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons") 

plot.by.site=ggplot(df.light, aes(Light)) + geom_density(aes(fill=Season), alpha=0.3, position = 'stack') +
  scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons") + facet_wrap(~Location, scales="free") + scale_fill_manual(values=c("darkorange1", "dodgerblue1"))

######
# can loop as a list
p=ggplot(df.light, aes(Light, fill=Season)) + geom_density(alpha=0.3, position = 'stack') + scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons") 

plots=dlply(df.light, .(Location), function(x) p %+% x + facet_wrap(~Location))


###### plot of light x depth by season with exponential curve fitting
Sum<-df.light[(df.light$Season=="summer"),]
Win<-df.light[(df.light$Season=="winter"),]

plot(Light~newDepth, Sum, col="coral", pch=16, xlab="Depth (m)", ylab="DLI at coral", 
     ylim=c(0, 30), xlim=c(0, 10))
summod<-lm(log(Light)~newDepth, Sum)
curve(exp(coef(summod)["(Intercept)"]+coef(summod)["newDepth"]*x), add=TRUE, col="coral", lwd=2)
par(new=T)
plot(Light~newDepth, Win, col="lightblue2", pch=16, xaxt="n", yaxt="n", 
     xlab="", ylab="", ylim=c(0, 30), xlim=c(0, 10))
wintmod<-lm(log(Light)~newDepth, Win)
curve(exp(coef(wintmod)["(Intercept)"]+coef(wintmod)["newDepth"]*x), add=TRUE, col="lightblue2", lwd=2)
legend("topright", legend=c("summer", "winter"), col=c("coral", "lightblue2"), lty=1, lwd=1, pch=16, bty="n")
**Figure** *Relationship between mean seasonal light availability (mean daily light integral) at depth of coral fragment collection*

Figure Relationship between mean seasonal light availability (mean daily light integral) at depth of coral fragment collection

Dissolved nutrients

Dissolved inorganic nutrients in seawater were quantified on water samples (ca. 100 ml) taken from each site and filtered through a GF/F filter and stored in acid-washed bottles and frozen at -20 °C until analyzes. Dissolved inorganic nutrients—ammonium (NH4+), nitrate + nitrite (NO3- + NO2-), phosphate (PO43-), and silicate (Si(OH)4)—were measured at University of Hawai‘i at Mānoa SOEST Laboratory for Analytical Biogeochemistry using a Seal Analytical AA3 HR nutrient autoanalyzer and expressed as μmol l-1.

Silicate showed no effects (p > 0.323). Phosphate increased in the winter (p=0.046) most notably at northern sites. N+N was higher in northern sites relative to southern sites and increased during the winter at all sites compared to the summer. Ammonium increased in the winter as well (p <0.001).

#########################
#########################
# nutrients
nutr<-read.csv("data/environmental/PanKBay_nutrients.csv")

#models
mod<-lm(silicate~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(phosphate~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(N.N~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(ammonium~Location+Date, data=nutr); Anova(mod, type=2)

nutr$Date<-mdy(nutr$Date) # corrects date format
dates<-cbind("10-Aug '16", "20-Dec '16")

RfHIMB<-nutr[(nutr$Location=="HIMB"), ]
Rf42<-nutr[(nutr$Location=="R42"), ]
Rf46<-nutr[(nutr$Location=="F1-46"), ]
RF10<-nutr[(nutr$Location=="F8-R10"), ]

Reefs<-c("Reef 46", "Reef 42", "Reef 10", "HIMB")
plot_colors<-c("mediumseagreen", "dodgerblue", "salmon", "orchid")

###############
# Phosphate
# par(mfrow=c(1,1), mar=c(2,4,1,1), mgp=c(2,0.5,0))

# figure layout
layout(matrix(c(1,2,3,4), nrow=1, byrow=TRUE))
par(mar=c(5,5,7,1))

###############
# Silicate
plot(y=Rf46$silicate, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("silicate"~(mu*mol~L^-1), sep="")), ylim=c(0, 15), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots HIMB
#plots Reef 42
with(Rf46, lines(Rf42$silicate, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots HIMB
with(RfHIMB, lines(RfHIMB$silicate, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots Reef 10
with(Rf46, lines(RF10$silicate, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))


###############
# Phosphate
plot(y=Rf46$phosphate, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("phosphate"~(mu*mol~L^-1), sep="")), ylim=c(0, 0.25), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots HIMB
#plots Reef 42
with(Rf46, lines(Rf42$phosphate, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots HIMB
with(RfHIMB, lines(RfHIMB$phosphate, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots Reef 10
with(Rf46, lines(RF10$phosphate, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))


###############
# N+N
plot(y=Rf46$N.N, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("nitrate+nitrite"~(mu*mol~L^-1), sep="")), ylim=c(0, 1.5), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots HIMB
#plots Reef 42
with(Rf46, lines(Rf42$N.N, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots HIMB
with(RfHIMB, lines(RfHIMB$N.N, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots Reef 10
with(Rf46, lines(RF10$N.N, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))


###############
# Ammonium
plot(y=Rf46$ammonium, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("ammonium"~(mu*mol~L^-1), sep="")), ylim=c(0, 1.5), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots HIMB
#plots Reef 42
with(Rf46, lines(Rf42$ammonium, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots HIMB
with(RfHIMB, lines(RfHIMB$ammonium, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots Reef 10
with(Rf46, lines(RF10$ammonium, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
legend("topleft", legend=Reefs, col=plot_colors[1:4], lwd=1, pch=19, lty=2, cex=0.8, bty="n", x.intersp=0.2, y.intersp=1.5, inset=c(0, 0), seg.len=1.5)

dev.copy(pdf, "figures/environmental/all.nutrients.pdf", encod="MacRoman", height=3, width=7)
dev.off()

Suspended particulates

Gravimetric SPM

Total (organic + inorganic) suspended particulate matter (SPM) was quantified by collecting and filtering 0.5 L of seawater from each location during each season. Seawater was sieved at 243 μm to remove large debris and filtered onto a pre-combusted (4h, 450°C) GF/F filter (0.7 μm). Salts were removed by light rinsing with ddH2O. Filters were dried at 60°C for 24h and re-weighed to closest 0.001 g to obtain total SPM. Masses were normalized to a liter of water.

###### SPM
spm<-read.csv("data/environmental/PanKBay_SPM.csv")
spm$SPM..g.l<-spm$SPM..g.l*1000 # change to mg
spm$Date<-factor(spm$Date, levels=c("8/10/16", "12/20/16"))
colnames(spm)[3]="SPM.mg"
mod<-lm(SPM.mg~Location+Date, data=spm); Anova(mod, type=2)
## Anova Table (Type II tests)
## 
## Response: SPM.mg
##           Sum Sq Df F value    Pr(>F)    
## Location  8.8382  3 10.3062 0.0003025 ***
## Date      1.4114  1  4.9373 0.0386252 *  
## Residuals 5.4312 19                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(mod), par.strip.text=list(cex=0.7))

rm(mean) # this will remove 'mean' from library and let the function below proceed.
mean<-aggregate(SPM.mg~Date + Location, spm, mean)
SE<-aggregate(SPM.mg~Date+Location, spm, std.error)

spm.df<-cbind(mean, SE[c(0,3)])
colnames(spm.df)=c("Date", "Location", "SPM", "SE")
spm.df$Date<-ordered(spm.df$Date, c("8/10/16", "12/20/16"))
spm.df$Location<-ordered(spm.df$Location, c("F1-46", "R42", "F8-10", "HIMB"))

Reefs<-c("Reef 46", "Reef 42", "Reef 10", "HIMB")
Date<-c("10-Aug '16", "20-Dec '16")
plot_colors2<-c("mediumturquoise", "skyblue3", "lightcoral", "plum")

#####
# figure formatting script
Fig.formatting<-(theme_classic()) +
  theme(text=element_text(size=10),
        axis.line=element_blank(),
        legend.justification=c(1,1), legend.position=c(1,1),
        legend.background = element_rect("transparent", colour=NA),
        legend.text=element_text(size=10),
        legend.title = element_blank(),
        panel.border = element_rect(fill=NA, colour = "black", size=0.8),
        aspect.ratio=1, 
        axis.ticks = element_line(colour = 'black', size = 0.4),
        axis.ticks.length=unit(0.25, "cm"),
        axis.text.y=element_text(
          margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10), 
        axis.text.x=element_text(
          margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10)) +
theme(legend.spacing.x = unit(0.2, 'cm'),
      legend.key.size = unit(0.4, "cm"))

pd <- position_dodge(0.7) #offset for error bars and columns
#####

# SPM (mg/l)
Fig.SPM<-ggplot(spm.df, aes(x=Date, y=SPM, fill=Location)) +
  geom_bar(colour=c("plum","lightcoral", "skyblue3","mediumturquoise","plum","lightcoral", "skyblue3", "mediumturquoise"),
           stat="identity", position = pd, width=0.7, size=0.4) +
  geom_errorbar(aes(ymin=SPM-SE, ymax=SPM+SE), 
                colour=c("plum","lightcoral", "skyblue3","mediumturquoise","plum","lightcoral", "skyblue3","mediumturquoise"),
                size=0.4, width=0, position=pd) +
  xlab("Sampling points") +
  ylab(expression(paste("SPM", ~(mg~L^-1), sep=""))) +
  scale_x_discrete(labels=Date) +
  scale_y_continuous(expand=c(0,0), limits=c(0, 5)) +
  scale_fill_manual(values=alpha(c(plot_colors2), 0.5), 
                    labels=Reefs)+ Fig.formatting; Fig.SPM

ggsave(file="figures/environmental/SPM.pdf", height=4, width=4)

Isotopes of SPM/plankton

Plankton tows and seawater collections to parameterized particulates and plankters as potential heterotrophic end members available to reef corals. Sampling was done at 4 sites where corals were collected [North: (Reef 42, fringe-Reef 46), and South: (HIMB, fringe-Reef 10)], as well as two reefs in central region of the bay where corals were not collected (Central: Reef 25, fringe-Reef 25) to increase spatial resolution of suspended particulate isotope sample sizes. Using size-fractioned materials collected in seawater and plankton tows, we examined the spatial (among reef sites) and temporal patterns (among seasons) in stable isotope values of size-fractioned end members. Carbon and nitrogen isotope values among the 6 reef sites were not significantly different (p > 0.140), season had marginal effects (p > 0.049), whereas fraction influenced isotope values significantly (p< 0.001). Therefore, isotope values were pooled among reefs and seasons to best interpret size-fraction isotope values. Samples were not acidified prior to analysis as to not affect nitrogen isotope values.

SWiso<-read.csv("data/isotopes_SW_all times.csv")

mod<-(lm(d13C~Location+Time.point+SW.fraction..um, data=SWiso)); Anova(mod, type=2) # no site effect
## Anova Table (Type II tests)
## 
## Response: d13C
##                  Sum Sq Df F value    Pr(>F)    
## Location         17.926  5  1.3418 0.2626013    
## Time.point        7.921  1  2.9645 0.0914188 .  
## SW.fraction..um  76.419  4  7.1504 0.0001286 ***
## Residuals       130.920 49                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod<-(lm(d15N~Location+Time.point+SW.fraction..um, data=SWiso)); Anova(mod, type=2) # no site effect
## Anova Table (Type II tests)
## 
## Response: d15N
##                  Sum Sq Df F value    Pr(>F)    
## Location         2.3255  5  1.7291   0.14557    
## Time.point       1.0935  1  4.0653   0.04927 *  
## SW.fraction..um 30.3773  4 28.2335 3.454e-12 ***
## Residuals       13.1802 49                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
SWiso$Location<-factor(SWiso$Location, levels=c("F1-46", "R42", "F2-R25", "R25", "F8-R10", "HIMB"))
SWiso$SW.fraction..um<-factor(SWiso$SW.fraction..um, levels=c(">243", "<243", "100-243", "10-100", "0-10"))

winter.data<-SWiso[(SWiso$Time.point=="winter"),]
summer.data<-SWiso[(SWiso$Time.point=="summer"),]

These box plots show the seasonal isotope values for carbon and nitrogen according to size fractions.

# plots of SW isotopes by Season-- first showing by size, then by site
op<-par(mfrow = c(2,2), mar=c(5,5,2,1))

######### Sizes
# Summer size d13C
plot(d13C~SW.fraction..um, data=summer.data, ylim=c(-25,-10), 
     main=expression(paste("summer"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
     col="paleturquoise3", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))

# Winter size d13C
plot(d13C~SW.fraction..um, data=winter.data, ylim=c(-25,-10), 
     main=expression(paste("winter"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
     col="paleturquoise3", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))

# Summer size d15N
plot(d15N~SW.fraction..um, data=summer.data, ylim=c(3,10), 
     main=expression(paste("summer"~ delta^{15}, "N")), col="plum", cex.axis=0.7, cex.main=1, cex.lab= 0.8,
     xlab=expression(paste("Size Fraction (", mu,m,")")),
     ylab=expression(paste(delta^{15}, "N (‰, air)")))

# Winter size d15N
plot(d15N~SW.fraction..um, data=winter.data, ylim=c(3,10), 
     main=expression(paste("winter"~ delta^{15}, "N")), ylab=expression(paste(delta^{15}, "N (‰, air)")),
     col="plum", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))
**Figure** *Size-fractioned sample isotope values stratified by seasons (summer and winter)*

Figure Size-fractioned sample isotope values stratified by seasons (summer and winter)

These figures show the same size-fractioned isotope data as above, but stratified by location and pooled across seasons.

######## Sites
op<-par(mfrow = c(2,2), mar=c(5,5,2,1))

# Summer site d13C
plot(d13C~Location, data=summer.data, ylim=c(-25,-10), 
     main=expression(paste("summer"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
     col="lightskyblue", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab="Reef Sites")

# Winter site d13C
plot(d13C~Location, data=winter.data, ylim=c(-25,-10), 
     main=expression(paste("winter"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
     col="lightskyblue", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab="Reef Sites")

# Summer site d15N
plot(d15N~Location, data=summer.data, ylim=c(3,10), 
     main=expression(paste("summer"~ delta^{15}, "N")), col="salmon", cex.axis=0.7, cex.main=1, cex.lab= 0.8,
     ylab=expression(paste(delta^{15}, "N (‰, air)")), xlab="Reef Sites")

# Winter site d15N
plot(d15N~Location, data=winter.data, ylim=c(3,10), 
     main=expression(paste("winter"~ delta^{15}, "N")), ylab=expression(paste(delta^{15}, "N (‰, air)")),
     col="salmon", cex.axis=0.7, cex.main=1, xlab="Reef Sites")
**Figure** *Isotope values in seawater particles stratified by sites (northwest to southeast)*

Figure Isotope values in seawater particles stratified by sites (northwest to southeast)

Pooling the isotope data across sites and seasons (where effects were absent or negligible), this planktonic foodweb for carbon and nitrogen illustrates differences in the food sources and their isotopic composition.

rm(mean)
####
#### making scatter for d15N and d13C, pooled across seasons and sites
mix.N.mean<-aggregate(d15N~SW.fraction..um, data=SWiso, mean)
mix.N.SE<-aggregate(d15N~SW.fraction..um, data=SWiso, std.error)
mix.C.mean<-aggregate(d13C~SW.fraction..um, data=SWiso, mean)
mix.C.SE<-aggregate(d13C~SW.fraction..um, data=SWiso, std.error)
mix.data<-cbind(mix.N.mean, mix.C.mean[c(2,0)], mix.N.SE[c(2,0)], mix.C.SE[c(2,0)]); colnames(mix.data)=c("fraction", "d15N", "d13C", "d15N.SE", "d13C.SE")

colors=c("#FF6A6A", "#00B2EE", "#FFB90F", "#3CB371", "#8B7500")
op<-par(mfrow = c(1,1), mar=c(5,4,1,5),xpd=TRUE, pty="sq")

size.labels=c(expression(paste("> 243"~mu,m)), expression(paste("< 243"~mu,m)), expression(paste("100-243"~mu,m)), expression(paste("10-100"~mu,m)), expression(paste("0-10"~mu,m)))

#### make the plot
plot(d15N~d13C, data=mix.data, type="n", ylim=c(4.5,8), xlim=c(-23,-17), tck=-0.03, cex.axis=0.7, cex.lab=0.7,
     ylab=expression(paste(delta^{15}, "N (‰, air)")), xlab=expression(paste(delta^{13}, "C (‰, V-PDB)")))
legend("topleft", inset=c(0.05,0.0), legend=size.labels, col=colors, pch=19, cex=0.6, bty="n", x.intersp=1.8, y.intersp = 1.4)
arrows(mix.data$d13C-mix.data$d13C.SE, mix.data$d15N, mix.data$d13C+mix.data$d13C.SE, mix.data$d15N, length=0.0, angle=90, code=3, lwd=0.7)
arrows(mix.data$d13C, mix.data$d15N-mix.data$d15N.SE, mix.data$d13C, mix.data$d15N+mix.data$d15N.SE, length=0.0, angle=90, code=3, lwd=0.7)
points(d15N~d13C, data=mix.data, pch=19, cex=0.8, col=colors)
**Figure** *Isotope values for size-fractioned seawater particles and plankters pooled across seasons and reef sites*

Figure Isotope values for size-fractioned seawater particles and plankters pooled across seasons and reef sites

dev.copy(pdf, "figures/environmental/iso.sources.KBay.pdf", encod="MacRoman", height=4, width=4)
dev.off()

Biology

################################
### Biological responses
################################

#data : this is the master file

# add in light at depth column from df.light dataframe
data$Light<-df.light$Light

##### produce a categorical depth bin ####
depth<-data$newDepth
data$depth.bin<-factor(ifelse(depth<2, "<2m", ifelse(depth >2 & depth <4, "2-4m", ifelse(depth >4 & depth <6, "4-6m", ">6m"))), levels=c("<2m", "2-4m", "4-6m", ">6m"))

aggregate(Sample.ID~depth.bin+Season+Location, data, length)
data$depth.bin.small<-factor(ifelse(depth<4, "<4m", ">4m"), levels= c("<4m", ">4m"))

################################################
# calculate, normalized dependent variables
################################################
str(data)
data$cells.ml<-as.numeric(data$cells.ml)

# helpful shorthand
SA<-data$surface.area # surface area in cm2
blastate<-data$total.blastate.ml # tissue slurry blastate in ml

# AFDW.mg. == convert AFDW g to mg, mutiply by blastate volume, divide by cm2
data$biomass<- (data$mg.biomass.ml*blastate)/SA

# Symbiodinium.cells. == cell.ml * blastate / SA
data$zoox<- (data$cells.ml*blastate)/SA

# total chlorophyll == ug.chl.a.ml * blastate + ug.chl.c2.ml * blastate / SA
data$chltot<-(data$ug.chl.a.ml)+(data$ug.chl.c2.ml)*blastate/SA

# pg.chlorophyll.a..cell + pg.chlorophyll.c2..cell == ug.chltot.ml * 10^6 / cells.ml
data$chlcell<- (data$ug.chl.a.ml*10^6+data$ug.chl.c2.ml*10^6)/data$cells.ml

qPCR

####################
# qPCR
#########

# qPCR
# Use steponeR to import data and calculate proporation of C and D symbionts
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path="data/qPCR", pattern = "txt$", full.names = T); Mcap.plates
Mcap <- steponeR(files=Mcap.plates, delim="\t",
                 target.ratios=c("C.D"),
                 fluor.norm=list(C=2.26827, D=0),
                 copy.number=list(C=33, D=3),
                 ploidy=list(C=1, D=1), 
                 extract=list(C=0.813, D=0.813))

Mcap <- Mcap$result
head(Mcap)

# remove +/-control
Mcap <- Mcap[grep("+C52", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("H2O", Mcap$Sample.Name, fixed=T, invert = T), ]

# to remove any early-amplification CT noise
Mcap$C.CT.mean[which(Mcap$C.CT.mean < 15)] <- 0

#Remove failed samples, i.e., those where either C or D were NOT found in both reps
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]

# replace CT means with 'NA' as zero
Mcap$C.CT.mean[is.na(Mcap$C.CT.mean)] <-0
Mcap$D.CT.mean[is.na(Mcap$D.CT.mean)] <-0

Mcap$C.D[is.na(Mcap$C.D)] <- 1 # sets all infinity (= 100% C) to 1.0

# caluclate proportion C and proprtion D where C and D are both present
Mcap$propC<- Mcap$C.D / (Mcap$C.D + 1)
Mcap$propD<- 1 / (Mcap$C.D + 1)

# where C and D are not cooccuring...
# if C.D = 1 = 100% C, make 'PropC' = 1 and 'PropD' = 0
# if C.D = 0 = 100% D, make 'PropD' = 1 and 'PropC' = 0
Mcap$propC[which(Mcap$C.D==1)] <- 1
Mcap$propD[which(Mcap$propC==1)] <- 0
Mcap$propD[which(Mcap$C.D==0)] <- 1

# calculate FOUR COMMUNITY categories: C, C>D, D>C, D
Mcap$Mix <- factor(ifelse(Mcap$propC > Mcap$propD, ifelse(Mcap$propD!= 0, "CD", "C"), ifelse(Mcap$propD > Mcap$propC, ifelse(Mcap$propC!=0, "DC", "D"), NA)), levels=c("C", "CD", "DC", "D"))

# Identify SINGLE dominant symbiont clade: C or D
Mcap$Dom <- factor(substr(as.character(Mcap$Mix), 1, 1))

# Set zeros to NA to facilitate log transformation
Mcap$propC[which(Mcap$propC==0)] <- NA
Mcap$propD[which(Mcap$propD==0)] <- NA

######## look for duplicates in dataset by year and type of event (bleach/recover)
Mcap[duplicated(Mcap$Sample.Name), ] ## duplicates

# remove duplicated
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_15" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_05" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R42_06" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_01" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_02" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_03" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_13" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_14" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_15" & Mcap$File.Name=="Wall_PanKbay_plate3.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_15" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),] 
Mcap<-Mcap[!(Mcap$Sample.Name=="R42_11" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),] 

# parse Sample ID and Site from "Site.Name"
Mcap<-cbind(Mcap, colsplit(Mcap$Sample.Name, pattern= "_", c("Location", "Sample.ID")))
Mcap$Season<-as.factor("winter")
Mcap$Location<-as.factor(Mcap$Location)

Mcap$Location<-revalue(Mcap$Location, c("R10"="F8-R10", "R46"="F1-R46")) # rename factor levels

# make new factors for bay region and reef type
Mcap$Bay.region <- ifelse(Mcap$Location=="R42" | Mcap$Location=="F1-R46", "northern", "southern")
Mcap$Reef.type <- ifelse(Mcap$Location=="R42" | Mcap$Location=="HIMB", "patch", "fringe")

### reorder columns and finish
Mcap<-Mcap[, c(17,15,19,18,16, 1:14)] # reordered to match masterdata, and finish

### structure winter and summer qPCR dataframes to have same columns and combine dataframe 
qPCR.winter<-Mcap[ , (names(Mcap) %in% 
            c("Season", "Location", "Reef.type", "Bay.region", "Sample.ID", "propC", "propD", "Mix", "Dom"))]
qPCR.summer<-qPCR.Innis[ , (names(qPCR.Innis) %in% 
            c("Season", "Location", "Reef.type", "Bay.region", "Sample.ID", "propC", "propD", "Mix", "Dom"))] 

# merge qPCR files
qPCR.all<-rbind(qPCR.winter, qPCR.summer)

# add to master data
data.all<-merge(data, qPCR.all, by=c("Season", "Location", "Reef.type", "Bay.region", "Sample.ID"), all.x=T)

###### remove columns no longer needed, update "Depth" to be tide-corrected depth )= newDepth
data.trim<-data.all[ , !(names(data.all) %in% c("total.blastate.ml", "Date", "Time.of.collection", "Depth..m", "Time.r", "surface.area.cm2", "cells.ml", "ug.chl.a.ml", "ug.chl.c2.ml", "mg.biomass.ml", "host..mass.mg", "host..ugN", "host..ugC", "symb..mass.mg", "symb..ugN", "symb..ugC", "stationId", "datum", "TimeUTC", "TideHT", "Time"))]

data.trim$symb..C.N[data.trim$symb..C.N>=12.520270]=NA # set this outlier to NA

qPCR figure The distribution of C and D symbiont types is modeled using a logistic regression with dominant symbiont community (C vs. D) as a response and Light-at-depth, Season, and reef Location as a predictor. A generalized linear model with a binomial distribution and logit link function is used. Model selection (by model AIC values) between a fully crossed model (crossed) and additive model (add) were compared. The additive mode best fit the data; main effects were tested in anova with a Chi-square test.

## qPCR figures

#####
model.data<-data.trim[,c(16, 1:4, 17:24, 6:15, 25:28)]

#########################
#########################
#########################

# Plots with Light at depth as the predictor
# qPCR and symbionts

#########################
###### Plot Dominant Symbiont and Depth (both seasons)
Symb<-model.data
Symb$Location <- factor(Symb$Location, levels = c("F1-R46", "R42", "F8-R10", "HIMB")) # set levels

Symb$Dominant <- ifelse(Symb$Dom=="D", 1, 0)
add <-glm(Dominant~Light+Season+Location, family = "binomial", data = Symb)
crossed <-glm(Dominant~Light*Season*Location, family = "binomial", data = Symb)
# anova(add, crossed) # additive model best
anova(add, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Dominant
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       118     140.51              
## Light     1  14.3551       117     126.16 0.0001514 ***
## Season    1   7.1646       116     118.99 0.0074356 ** 
## Location  3   8.2623       113     110.73 0.0408905 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
results.all<-add

#summary(results.all)
plot(allEffects(results.all))
posthoc<-emmeans(results.all, ~Location, type = "response") # adding 'response' backtransforms data
CLD(posthoc, Letters=letters)
##  Location       prob         SE  df  asymp.LCL asymp.UCL .group
##  R42      0.07342293 0.04408130 Inf 0.02176986 0.2200623  a    
##  F1-R46   0.19712898 0.07197537 Inf 0.09148303 0.3744878  a    
##  F8-R10   0.34040087 0.11812381 Inf 0.15542876 0.5913683  a    
##  HIMB     0.34225969 0.09850311 Inf 0.18081339 0.5509143  a    
## 
## Results are averaged over the levels of: Season 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale 
## significance level used: alpha = 0.05
results.all.mod <-glm(Dominant~Light, family = "binomial", data = Symb)
fitted.all <- predict(results.all.mod, newdata = list(Light=seq(0,25,0.1)), type = "response")

###### summer only symbionts
sum.dat<-Symb[(Symb$Season=="summer"),]
results.sum=glm(Dominant~Light, family = "binomial", data = sum.dat)
anova(results.sum, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Dominant
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                     59      67.48              
## Light  1    16.95        58      50.53 3.838e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(results.sum)
fitted.sum <- predict(results.sum, newdata = list(Light=seq(0,25,0.1)), type = "response")
# plot(fitted.sum, ylab="proportion D", ylim=c(0,1))

###### winter only symbionts
wint.dat<-Symb[(Symb$Season=="winter"),]
results.win=glm(Dominant~Light, family = "binomial", data = wint.dat)
anova(results.win, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Dominant
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                     58     72.583           
## Light  1   4.8011        57     67.782  0.02844 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(results.win)
fitted.win <- predict(results.win, newdata = list(Light=seq(0,25,0.1)), type = "response")
# plot(fitted.win, ylab="proportion D", ylim=c(0,1))

Logistic regression of the probability of Montipora capitata have a clade D-dominant symbiont community (probability of 1.0) versus one dominated by clade C symbiont (probability of 0). Both seasons show similar patterns, although the probability of corals being clade D-dominant increases in winter months as a function of light-at-depth relative to summer months. This likely is explained by much less light in winter months and a reduced magnitude in the attenuation gradient across depth.

######
######
## Figure of Dominant symbiont clades across seasons
## **Note** where points equal 0.0 is 100% D, where they equal 0 is 100% C
par(mar=c(5,4,3,2))
plot(sum.dat$propD~sum.dat$Light, xlab="Light (DLI)", ylab = "Proportion of Clade D Symbiont", pch=19, col="salmon", xlim=c(0,25), ylim=c(0.0, 1.0), cex=0.5)
par(new=T)
plot(wint.dat$propD~wint.dat$Light, pch=19, col="lightskyblue", xlim=c(0,25), ylim=c(0.0, 1.0), xaxt="n", xlab="", ylab="", cex=0.5)
lines(fitted.all ~ seq(0,25,0.1), col="gray30", lwd=1, lty=2)
lines(fitted.sum ~ seq(0,25,0.1), col="coral", lwd=1)
lines(fitted.win ~ seq(0,25,0.1), col="lightskyblue", lwd=1)
legend("bottomright", pch=c(19,19, NA), lty=c(1,1,2), col=c("coral", "lightskyblue", "gray30"), legend=c("Summer", "Winter", "Combined"), lwd=1, bty="n", x.intersp=1, pt.cex=0.8, y.intersp=2, cex=0.7, inset=c(0.0, 0.15), seg.len=3)
**Figure** Probability of *Montipora capitata* hosting clade D symbionts across a light availability (DLI) gradient in summer and winter seasons. Lines represent model fit to data by logistic regression for summer alone, winter alone, and combined (summer + winter) data

Figure Probability of Montipora capitata hosting clade D symbionts across a light availability (DLI) gradient in summer and winter seasons. Lines represent model fit to data by logistic regression for summer alone, winter alone, and combined (summer + winter) data

dev.copy(pdf, "figures/symbionts/Symbionts_by_Season_light.pdf", encod="MacRoman", height=4, width=4)
dev.off()

Seen another way–Histograms. By separating these data into 2 histograms the probability of clade D- (1.0) and clade C-dominance (0.0) can be visualized. Here it is clear that there is a greater probability of finding corals with clade D dominant communites in the shallows during summer times, whereas in winter there is a slightly greater probability of finding clade D in corals with increasing depth/low light.

#######
#######
par(mfrow=c(1,2))
Dom1 <- subset(Symb, !is.na(newDepth) & !is.na(Dominant))
Dom.sum<-Dom1[(Dom1$Season=="summer"),]
Dom.win<-Dom1[(Dom1$Season=="winter"),]

logi.hist.plot(Dom.sum$Light, Dom.sum$Dominant, boxp = FALSE, type = "hist", col="coral", xlabel = "Light (DLI)", ylabel = "", ylabel2 = "", main="summer")
mtext(side = 2, text = "Probability of Clade D", line = 3, cex = 1)

logi.hist.plot(Dom.win$Light, Dom.win$Dominant, boxp = FALSE, type = "hist", col="lightskyblue", xlabel = "Light (DLI)", ylabel = "", ylabel2 = "", main="winter")
mtext(side = 4, text = "Frequency", line = 0.5, cex=1)
mtext(side = 4, text = "C                                             D", line = 0.5, cex = 0.8)
**Figure** Histograms of probability of *Montipora capitata* hosting clade D symbionts across a gradient in light availability (DLI) among seasons. Lines represent model fit to data by logistic regression

Figure Histograms of probability of Montipora capitata hosting clade D symbionts across a gradient in light availability (DLI) among seasons. Lines represent model fit to data by logistic regression

dev.copy(pdf, "figures/symbionts/Symb_Season_Light_logistic.pdf", encod="MacRoman", height=5, width=8)
dev.off()

If we use Depth (m) as a predictor, main effects in logstic regression (i.e., glm) are comparable to logistic regressions with respect to light-at-depth (effect of Depth in glm p <0.001), but seasonal effects are absent (Season effect is p = 0.747).

Predictions of depth effects on the proability of clade D dominance using a glm with depth alone in both seasons, together and separately, are seen in the figure.

#########################
#########################

# this code runs the same as above, except using colony depth as a predictor instead of light-at-depth


# Depth as a predictor
# qPCR and symbionts

#########################
###### Plot Dominant Symbiont and Depth (both seasons)
Symb<-model.data
Symb$Dominant <- ifelse(Symb$Dom=="D", 1, 0)
add<-glm(Dominant~newDepth+Season+Location, family = "binomial", data = Symb)
crossed<-glm(Dominant~newDepth*Season*Location, family = "binomial", data = Symb)
#AIC(add, crossed)

anova(add, test = "Chisq")
#summary(add)

results.all<-add

results.all.mod<-glm(Dominant~newDepth, family = "binomial", data = Symb)
fitted.all <- predict(results.all.mod, newdata = list(newDepth=seq(0,10,0.1)), type = "response")

###### summer only symbionts
sum.dat<-Symb[(Symb$Season=="summer"),]
results.sum<-glm(Dominant~newDepth, family = "binomial", data = sum.dat)
#anova(results.sum, test = "Chisq")
#summary(results.sum)
fitted.sum <- predict(results.sum, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
# plot(fitted.sum, ylab="proportion D", ylim=c(0,1))

###### winter only symbionts
wint.dat<-Symb[(Symb$Season=="winter"),]
results.win<-glm(Dominant~newDepth, family = "binomial", data = wint.dat)
#anova(results.win, test = "Chisq")
#summary(results.win)
fitted.win <- predict(results.win, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
# plot(fitted.win, ylab="proportion D", ylim=c(0,1))

######
######
## Figure of Dominant symbiont clades across seasons
## **Note** where points equal 0.0 is 100% D, where they equal 0 is 100% C
par(mar=c(5,4,3,2))
plot(sum.dat$propD~sum.dat$newDepth, xlab="Depth (m)", ylab = "Proportion of Clade D Symbiont", pch=19, col="salmon", xlim=c(0,10), ylim=c(0.0, 1.0), cex=0.5)
par(new=T)
plot(wint.dat$propD~wint.dat$newDepth, pch=19, col="lightskyblue", xlim=c(0,10), ylim=c(0.0, 1.0), xaxt="n", xlab="", ylab="", cex=0.5)
lines(fitted.all ~ seq(0,10,0.1), col="gray30", lwd=1, lty=2)
lines(fitted.sum ~ seq(0,10,0.1), col="coral", lwd=1)
lines(fitted.win ~ seq(0,10,0.1), col="lightskyblue", lwd=1)
legend("topright", pch=c(19,19, NA), lty=c(1,1,2), col=c("coral", "lightskyblue", "gray30"), legend=c("Summer", "Winter", "Combined"), lwd=1, bty="n", x.intersp=1, pt.cex=0.7, y.intersp=2, cex=0.5, inset=c(0, 0.15), seg.len=3)
**Figure** Probability of *Montipora capitata* hosting clade D symbionts across a **DEPTH gradient** in summer and winter seasons. Lines represent model fit to data by logistic regression for summer alone, winter alone, and combined (summer + winter) data

Figure Probability of Montipora capitata hosting clade D symbionts across a DEPTH gradient in summer and winter seasons. Lines represent model fit to data by logistic regression for summer alone, winter alone, and combined (summer + winter) data

dev.copy(pdf, "figures/symbionts/Symbionts_by_Season.pdf", encod="MacRoman", height=4, width=4)
dev.off()
#######
#######
par(mfrow=c(1,2))
Dom1 <- subset(Symb, !is.na(newDepth) & !is.na(Dominant))
Dom.sum<-Dom1[(Dom1$Season=="summer"),]
Dom.win<-Dom1[(Dom1$Season=="winter"),]

logi.hist.plot(Dom.sum$newDepth, Dom.sum$Dominant, boxp = FALSE, type = "hist", col="coral", xlabel = "Depth (m)", ylabel = "", ylabel2 = "", main="summer")
mtext(side = 2, text = "Probability of Clade D", line = 3, cex = 1)

logi.hist.plot(Dom.win$newDepth, Dom.win$Dominant, boxp = FALSE, type = "hist", col="lightskyblue", xlabel = "Depth (m)", ylabel = "", ylabel2 = "", main="winter")
mtext(side = 4, text = "Frequency", line = 0.5, cex=1)
mtext(side = 4, text = "C                                             D", line = 0.5, cex = 0.8)

dev.copy(pdf, "figures/symbionts/Symb_Season_Depth_logistic.pdf", encod="MacRoman", height=5, width=8)
dev.off()

Model assumptions

Physiology models

Total Biomass

########## ########## 
######### biomass ---- 
Y<-model.data$biomass
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
summary(full)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season * Light * Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 811
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.81107 -0.71811 -0.02691  0.62021  3.10297 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept)  7.214   2.686   
##  Residual             59.912   7.740   
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##                           Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)              2.742e+01  2.581e+00  2.116e+01  10.624  6.1e-10
## Seasonwinter            -5.409e+00  3.111e+00  1.088e+02  -1.738    0.085
## Light                   -1.969e-01  2.245e-01  1.107e+02  -0.877    0.382
## DomD                     1.722e+00  8.125e+00  1.095e+02   0.212    0.833
## Seasonwinter:Light       7.319e-01  4.225e-01  1.106e+02   1.732    0.086
## Seasonwinter:DomD        2.354e+00  9.431e+00  1.095e+02   0.250    0.803
## Light:DomD               9.707e-04  5.473e-01  1.093e+02   0.002    0.999
## Seasonwinter:Light:DomD -6.400e-02  8.240e-01  1.097e+02  -0.078    0.938
##                            
## (Intercept)             ***
## Seasonwinter            .  
## Light                      
## DomD                       
## Seasonwinter:Light      .  
## Seasonwinter:DomD          
## Light:DomD                 
## Seasonwinter:Light:DomD    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD   Ssnw:L Ssn:DD Lgh:DD
## Seasonwintr -0.566                                          
## Light       -0.727  0.550                                   
## DomD        -0.214  0.190  0.207                            
## Ssnwntr:Lgh  0.330 -0.818 -0.454 -0.129                     
## Ssnwntr:DmD  0.182 -0.333 -0.175 -0.865  0.277              
## Light:DomD   0.281 -0.236 -0.386 -0.943  0.211  0.815       
## Ssnwnt:L:DD -0.172  0.416  0.238  0.636 -0.510 -0.873 -0.671
print(anova(full, type=2), digits=4)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##                  Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)  
## Season           184.52  184.52     1 108.7   3.080 0.0821 .
## Light              2.54    2.54     1 105.8   0.042 0.8372  
## Dom               43.22   43.22     1 108.3   0.721 0.3975  
## Season:Light     232.07  232.07     1 110.5   3.873 0.0516 .
## Season:Dom         3.73    3.73     1 109.5   0.062 0.8034  
## Light:Dom          0.27    0.27     1 108.8   0.005 0.9463  
## Season:Light:Dom   0.36    0.36     1 109.7   0.006 0.9382  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ranef(full)
## $Location
##        (Intercept)
## F1-R46    1.740698
## F8-R10    1.887255
## HIMB     -3.083176
## R42      -0.544777
fixef(full)
##             (Intercept)            Seasonwinter                   Light 
##             27.41887497             -5.40870737             -0.19690467 
##                    DomD      Seasonwinter:Light       Seasonwinter:DomD 
##              1.72219130              0.73193108              2.35376651 
##              Light:DomD Seasonwinter:Light:DomD 
##              0.00097068             -0.06400212
sjp.lmer(full, y.offset = .4)

#sjp.lmer(full, vars = "Seasonwinter", type = "ri.slope")

posthoc<-emmeans(full, ~Light:Season)
CLD(posthoc, Letters=letters)
##     Light Season   emmean       SE    df lower.CL upper.CL .group
##  8.043006 summer 26.70017 2.518114 18.90 21.42787 31.97246  a    
##  8.043006 winter 28.09789 1.822194  5.64 23.56837 32.62741  a    
## 
## Results are averaged over the levels of: Dom 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="biomass", par.strip.text=list(cex=0.7))

Chlorophyll a

######### chltotal-- 
Y<-model.data$chltot
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
summary(full)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season * Light * Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 419.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.43730 -0.62654 -0.04514  0.46171  3.03177 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.3024   0.5499  
##  Residual             1.7520   1.3236  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)               5.50685    0.46704  15.02000  11.791 5.40e-09
## Seasonwinter              2.67199    0.53227 108.49000   5.020 2.04e-06
## Light                    -0.05151    0.03848 110.33000  -1.339   0.1834
## DomD                      0.80228    1.39071 109.08000   0.577   0.5652
## Seasonwinter:Light       -0.15910    0.07242 110.22000  -2.197   0.0301
## Seasonwinter:DomD        -3.36688    1.61430 109.06000  -2.086   0.0393
## Light:DomD               -0.10411    0.09366 108.91000  -1.111   0.2688
## Seasonwinter:Light:DomD   0.23738    0.14107 109.24000   1.683   0.0953
##                            
## (Intercept)             ***
## Seasonwinter            ***
## Light                      
## DomD                       
## Seasonwinter:Light      *  
## Seasonwinter:DomD       *  
## Light:DomD                 
## Seasonwinter:Light:DomD .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD   Ssnw:L Ssn:DD Lgh:DD
## Seasonwintr -0.533                                          
## Light       -0.689  0.547                                   
## DomD        -0.201  0.190  0.205                            
## Ssnwntr:Lgh  0.308 -0.818 -0.447 -0.130                     
## Ssnwntr:DmD  0.171 -0.333 -0.173 -0.865  0.277              
## Light:DomD   0.265 -0.235 -0.385 -0.943  0.210  0.815       
## Ssnwnt:L:DD -0.161  0.417  0.235  0.636 -0.510 -0.873 -0.671
print(anova(full, type=2), digits=4)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##                  Sum Sq Mean Sq NumDF DenDF F.value   Pr(>F)    
## Season            36.79   36.79     1 108.4  20.999 1.24e-05 ***
## Light             18.40   18.40     1 109.0  10.502  0.00158 ** 
## Dom               10.32   10.32     1 108.1   5.891  0.01687 *  
## Season:Light       4.34    4.34     1 110.1   2.475  0.11856    
## Season:Dom         7.62    7.62     1 109.1   4.350  0.03934 *  
## Light:Dom          0.00    0.00     1 108.5   0.000  0.98600    
## Season:Light:Dom   4.96    4.96     1 109.2   2.831  0.09530 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(full, ~Dom:Light:Season)
CLD(posthoc, Letters=letters)
##  Dom    Light Season   emmean        SE    df lower.CL upper.CL .group
##  D   8.043006 winter 4.992202 0.4220449 10.64 4.059429 5.924976  a    
##  D   8.043006 summer 5.057493 0.7525279 62.15 3.553282 6.561704  ab   
##  C   8.043006 summer 5.092535 0.3386960  4.61 4.199143 5.985928  a    
##  C   8.043006 winter 6.484888 0.3917497  7.68 5.574853 7.394923   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="total chlor", par.strip.text=list(cex=0.7))

ranef(full)
## $Location
##        (Intercept)
## F1-R46  -0.2528478
## F8-R10  -0.5465092
## HIMB     0.2274066
## R42      0.5719504
fixef(full)
##             (Intercept)            Seasonwinter                   Light 
##              5.50685274              2.67199375             -0.05151277 
##                    DomD      Seasonwinter:Light       Seasonwinter:DomD 
##              0.80227919             -0.15909984             -3.36688095 
##              Light:DomD Seasonwinter:Light:DomD 
##             -0.10410559              0.23737865
sjp.lmer(full, y.offset = .4)

Chlorophyll per symbiont cell

######### chlcell -- 
Y<-model.data$chlcell
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     6 4.8574 21.532 3.5713  -7.1426                         
## object 10 6.9000 34.691 6.5500 -13.1000 5.9574      4     0.2023
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 13.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7400 -0.5161 -0.0852  0.5185  2.4510 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.02551  0.1597  
##  Residual             0.05208  0.2282  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    1.87150    0.09786   5.69000  19.124 2.23e-06 ***
## Seasonwinter   0.07957    0.04766 112.84000   1.669   0.0978 .  
## Light         -0.03471    0.00525 114.62000  -6.611 1.26e-09 ***
## DomD          -0.52057    0.05268 112.94000  -9.881  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.407              
## Light       -0.478  0.476       
## DomD         0.106 -0.261 -0.447
print(anova(add, type=2), digits=5)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##        Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Season 0.1452  0.1452     1 112.84   2.787   0.09779 .  
## Light  2.2767  2.2767     1 114.62  43.711 1.261e-09 ***
## Dom    5.0854  5.0854     1 112.94  97.637 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(add), ylab="chla cell", par.strip.text=list(cex=0.7))

ranef(add)
## $Location
##         (Intercept)
## F1-R46  0.002272376
## F8-R10 -0.216768588
## HIMB    0.082748226
## R42     0.131747985
rand(add)
## Analysis of Random effects Table:
##          Chi.sq Chi.DF p.value    
## Location   23.9      1   1e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fixef(add)
##  (Intercept) Seasonwinter        Light         DomD 
##   1.87150224   0.07956797  -0.03471006  -0.52057321
sjp.lmer(full, y.offset = .4)

Isotope models

host d13C

########## host..d13C -- 
Y<-model.data$host..d13C
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     6 323.06 339.73 -155.53   311.06                         
## object 10 326.01 353.80 -153.00   306.01 5.0527      4     0.2819
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 321.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.85036 -0.66017 -0.07034  0.63948  2.03805 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.2418   0.4918  
##  Residual             0.7654   0.8749  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -16.37006    0.32741   7.43000 -49.998 1.14e-10 ***
## Seasonwinter   0.06996    0.18251 113.27000   0.383    0.702    
## Light          0.11577    0.02004 115.00000   5.777 6.62e-08 ***
## DomD          -1.12836    0.20172 113.40000  -5.594 1.56e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.465              
## Light       -0.546  0.474       
## DomD         0.119 -0.259 -0.445
print(anova(add, type=2), digits=5)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##         Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Season  0.1124  0.1124     1 113.27   0.147    0.7022    
## Light  25.5408 25.5408     1 115.00  33.369 6.615e-08 ***
## Dom    23.9479 23.9479     1 113.40  31.288 1.561e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(add), ylab="d13Chost", par.strip.text=list(cex=0.7))

ranef(add)
## $Location
##        (Intercept)
## F1-R46  -0.3582378
## F8-R10   0.4633780
## HIMB    -0.4407256
## R42      0.3355854
rand(add)
## Analysis of Random effects Table:
##          Chi.sq Chi.DF p.value    
## Location   17.9      1   2e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fixef(add)
##  (Intercept) Seasonwinter        Light         DomD 
## -16.37005799   0.06995578   0.11577346  -1.12836185
sjp.lmer(full, y.offset = .4)

symbiont d13C

########## symb..d13C --
Y<-model.data$symb..d13C
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## ..1     6 340.93 357.61 -164.47   328.93                           
## object 10 337.87 365.66 -158.94   317.87 11.061      4    0.02589 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 338.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.95192 -0.54508  0.00165  0.69202  2.68767 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.3278   0.5725  
##  Residual             0.8851   0.9408  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -16.59246    0.36895   6.73000 -44.972 1.34e-09 ***
## Seasonwinter  -0.02007    0.19635 113.11000  -0.102 0.918758    
## Light          0.12245    0.02159 114.92000   5.672 1.07e-07 ***
## DomD          -0.78673    0.21703 113.23000  -3.625 0.000434 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.444              
## Light       -0.522  0.475       
## DomD         0.115 -0.260 -0.446
print(anova(add, type=2), digits=5)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##         Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Season  0.0092  0.0092     1 113.11   0.010 0.9187576    
## Light  28.4743 28.4743     1 114.92  32.170  1.07e-07 ***
## Dom    11.6309 11.6309     1 113.23  13.141 0.0004344 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(add), ylab="d13Csymb", par.strip.text=list(cex=0.7))

ranef(add)
## $Location
##        (Intercept)
## F1-R46  -0.4229777
## F8-R10   0.6641010
## HIMB    -0.4709992
## R42      0.2298759
rand(add)
## Analysis of Random effects Table:
##          Chi.sq Chi.DF p.value    
## Location   20.7      1   5e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fixef(add)
##  (Intercept) Seasonwinter        Light         DomD 
## -16.59246334  -0.02007215   0.12244774  -0.78672936
sjp.lmer(full, y.offset = .4)

host-symbiont d13C

######### d13C..host.sym --
Y<-model.data$d13C..host.sym
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Light +Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Light + Season:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     8 28.734 50.967 -6.3670   12.734                         
## object 10 32.339 60.131 -6.1697   12.339 0.3946      2     0.8209
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: 
## Y ~ Season + Light + Dom + Season:Light + Season:Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 43
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0263 -0.4446  0.0338  0.5208  3.0008 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.01887  0.1374  
##  Residual             0.06362  0.2522  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)          0.133392   0.096724   8.690000   1.379  0.20233   
## Seasonwinter         0.296870   0.092051 110.200000   3.225  0.00166 **
## Light               -0.004318   0.006787 112.000000  -0.636  0.52598   
## DomD                -0.113913   0.087914 110.560000  -1.296  0.19776   
## Seasonwinter:Light  -0.024891   0.011636 111.280000  -2.139  0.03461 * 
## Seasonwinter:DomD   -0.288488   0.115782 110.470000  -2.492  0.01420 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD   Ssnw:L
## Seasonwintr -0.484                            
## Light       -0.587  0.554                     
## DomD         0.137 -0.111 -0.517              
## Ssnwntr:Lgh  0.284 -0.778 -0.484  0.252       
## Ssnwntr:DmD -0.086  0.036  0.364 -0.750 -0.357
print(anova(add, type=2), digits=4)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##              Sum Sq Mean Sq NumDF DenDF F.value   Pr(>F)    
## Season       0.7091  0.7091     1 110.2  11.147  0.00115 ** 
## Light        0.2799  0.2799     1 113.0   4.400  0.03817 *  
## Dom          1.4788  1.4788     1 111.4  23.246 4.55e-06 ***
## Season:Light 0.2911  0.2911     1 111.3   4.576  0.03461 *  
## Season:Dom   0.3950  0.3950     1 110.5   6.208  0.01420 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(add), ylab="d13C-hs", par.strip.text=list(cex=0.7))

ranef(add)
## $Location
##        (Intercept)
## F1-R46  0.04988743
## F8-R10 -0.19135411
## HIMB    0.04755924
## R42     0.09390744
rand(add)
## Analysis of Random effects Table:
##          Chi.sq Chi.DF p.value    
## Location   12.4      1   4e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fixef(add)
##        (Intercept)       Seasonwinter              Light 
##        0.133392463        0.296870058       -0.004317635 
##               DomD Seasonwinter:Light  Seasonwinter:DomD 
##       -0.113913174       -0.024890826       -0.288487575
sjp.lmer(add, y.offset = .4)

skeleton d13C

####### d13C..skel --
Y<-model.data$d13C..skel
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     6 311.33 328.01 -149.67   299.33                         
## object 10 315.61 343.40 -147.81   295.61 3.7238      4     0.4447
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 310.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96903 -0.81931  0.00245  0.67065  2.42406 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.1369   0.3700  
##  Residual             0.7036   0.8388  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   -2.839026   0.277165  10.130000 -10.243 1.15e-06 ***
## Seasonwinter   0.461515   0.174708 113.770000   2.642  0.00941 ** 
## Light          0.006456   0.019089 114.340000   0.338  0.73583    
## DomD           0.044872   0.193058 113.950000   0.232  0.81662    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.523              
## Light       -0.614  0.471       
## DomD         0.132 -0.257 -0.442
print(anova(add, type=2), digits=4)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##        Sum Sq Mean Sq NumDF DenDF F.value  Pr(>F)   
## Season  4.910   4.910     1 113.8   6.978 0.00941 **
## Light   0.080   0.080     1 114.3   0.114 0.73583   
## Dom     0.038   0.038     1 114.0   0.054 0.81662   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(add, ~Season)
CLD(posthoc, Letters=letters)
##  Season    emmean        SE   df  lower.CL  upper.CL .group
##  summer -2.764663 0.2270645 4.74 -3.358112 -2.171214  a    
##  winter -2.303148 0.2187764 4.14 -2.902642 -1.703655   b   
## 
## Results are averaged over the levels of: Dom 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="d13C-skel", par.strip.text=list(cex=0.7))

host d15N

####### host..d15N --
Y<-model.data$host..d15N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+Season:Light+ Season:Dom+ Light:Dom +(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Light + Season:Dom + Light:Dom + 
## ..1:     (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     9 86.794 111.81 -34.397   68.794                         
## object 10 88.753 116.55 -34.377   68.753 0.0406      1     0.8404
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: 
## Y ~ Season + Light + Dom + Season:Light + Season:Dom + Light:Dom +  
##     (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 100.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.85472 -0.56530  0.00305  0.63356  2.24611 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.29222  0.5406  
##  Residual             0.09533  0.3088  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          4.873694   0.283900   3.570000  17.167 0.000148 ***
## Seasonwinter         0.200690   0.113010 109.040000   1.776 0.078544 .  
## Light               -0.016754   0.008801 109.290000  -1.904 0.059587 .  
## DomD                -0.389112   0.250457 109.030000  -1.554 0.123176    
## Seasonwinter:Light  -0.022145   0.014615 109.190000  -1.515 0.132592    
## Seasonwinter:DomD    0.205938   0.183732 109.070000   1.121 0.264808    
## Light:DomD           0.032102   0.016224 109.040000   1.979 0.050375 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD   Ssnw:L Ssn:DD
## Seasonwintr -0.195                                   
## Light       -0.259  0.499                            
## DomD        -0.049 -0.107  0.073                     
## Ssnwntr:Lgh  0.096 -0.773 -0.373  0.293              
## Ssnwntr:DmD  0.024  0.069  0.065 -0.822 -0.402       
## Light:DomD   0.082  0.067 -0.315 -0.903 -0.210  0.635
print(anova(add, type=2), digits=4)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##              Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)  
## Season       0.2754  0.2754     1 109.0   2.888 0.0921 .
## Light        0.5038  0.5038     1 109.7   5.284 0.0234 *
## Dom          0.1139  0.1139     1 109.1   1.195 0.2767  
## Season:Light 0.2189  0.2189     1 109.2   2.296 0.1326  
## Season:Dom   0.1198  0.1198     1 109.1   1.256 0.2648  
## Light:Dom    0.3732  0.3732     1 109.0   3.915 0.0504 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(full), ylab="d15Nhost", par.strip.text=list(cex=0.7))

ranef(add)
## $Location
##        (Intercept)
## F1-R46  0.39986183
## F8-R10  0.03118106
## HIMB    0.33770089
## R42    -0.76874378
rand(add)
## Analysis of Random effects Table:
##          Chi.sq Chi.DF p.value    
## Location    120      1  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fixef(add)
##        (Intercept)       Seasonwinter              Light 
##         4.87369362         0.20069043        -0.01675445 
##               DomD Seasonwinter:Light  Seasonwinter:DomD 
##        -0.38911230        -0.02214507         0.20593819 
##         Light:DomD 
##         0.03210159
sjp.lmer(add, y.offset = .4)

symbiont d15N

######## symb..d15N --
Y<-model.data$symb..d15N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+Season:Light+ Season:Dom+ Light:Dom + (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Light + Season:Dom + Light:Dom + 
## ..1:     (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     9 97.046 122.06 -39.523   79.046                         
## object 10 98.398 126.19 -39.199   78.398 0.6483      1     0.4207
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: 
## Y ~ Season + Light + Dom + Season:Light + Season:Dom + Light:Dom +  
##     (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 110.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2292 -0.5979  0.1164  0.6161  1.9877 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.4345   0.6592  
##  Residual             0.1028   0.3206  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          4.344328   0.341713   3.410000  12.713 0.000546 ***
## Seasonwinter        -0.065757   0.117357 109.030000  -0.560 0.576416    
## Light               -0.004593   0.009141 109.210000  -0.502 0.616330    
## DomD                -0.390860   0.260091 109.020000  -1.503 0.135786    
## Seasonwinter:Light  -0.002555   0.015178 109.140000  -0.168 0.866608    
## Seasonwinter:DomD    0.488751   0.190804 109.050000   2.562 0.011787 *  
## Light:DomD           0.027628   0.016848 109.030000   1.640 0.103925    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD   Ssnw:L Ssn:DD
## Seasonwintr -0.169                                   
## Light       -0.224  0.499                            
## DomD        -0.042 -0.107  0.073                     
## Ssnwntr:Lgh  0.083 -0.773 -0.372  0.293              
## Ssnwntr:DmD  0.021  0.069  0.065 -0.822 -0.402       
## Light:DomD   0.071  0.067 -0.315 -0.903 -0.210  0.635
print(anova(add, type=2), digits=3)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##              Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)  
## Season        0.057   0.057     1   109    0.56  0.457  
## Light         0.001   0.001     1   110    0.01  0.937  
## Dom           0.124   0.124     1   109    1.21  0.274  
## Season:Light  0.003   0.003     1   109    0.03  0.867  
## Season:Dom    0.675   0.675     1   109    6.56  0.012 *
## Light:Dom     0.276   0.276     1   109    2.69  0.104  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(full, ~Season*Dom)
CLD(posthoc, Letters=letters)
##  Season Dom   emmean        SE   df lower.CL upper.CL .group
##  winter C   4.206022 0.3338188 3.20 3.180359 5.231685  a    
##  summer D   4.210158 0.3683973 4.74 3.247190 5.173125  ab   
##  summer C   4.306919 0.3304208 3.07 3.269372 5.344466  a    
##  winter D   4.546253 0.3360060 3.29 3.527659 5.564846   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="d15Nsymb", par.strip.text=list(cex=0.7))

ranef(add)
## $Location
##        (Intercept)
## F1-R46   0.4288524
## F8-R10   0.1246317
## HIMB     0.4085551
## R42     -0.9620392
rand(add)
## Analysis of Random effects Table:
##          Chi.sq Chi.DF p.value    
## Location    144      1  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fixef(add)
##        (Intercept)       Seasonwinter              Light 
##        4.344327831       -0.065756674       -0.004593471 
##               DomD Seasonwinter:Light  Seasonwinter:DomD 
##       -0.390860403       -0.002555484        0.488750689 
##         Light:DomD 
##        0.027627765
sjp.lmer(add, y.offset = .4)

host-symbiont d15N

####### d15N..host.sym  -- 
Y<-model.data$d15N..host.sym
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom+  (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     7 1.6903 21.144 6.1548  -12.310                         
## object 10 3.5687 31.360 8.2156  -16.431 4.1216      3     0.2486
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: 12
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.80715 -0.73247  0.06926  0.71888  2.22092 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Location (Intercept) 0.01318  0.1148  
##  Residual             0.05131  0.2265  
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         0.555024   0.080132   8.180000   6.926 0.000109 ***
## Seasonwinter        0.155035   0.051874 112.000000   2.989 0.003444 ** 
## Light              -0.015260   0.005323 113.890000  -2.867 0.004942 ** 
## DomD                0.090965   0.076384 112.290000   1.191 0.236207    
## Seasonwinter:DomD  -0.380440   0.097092 111.940000  -3.918 0.000154 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light  DomD  
## Seasonwintr -0.454                     
## Light       -0.556  0.322              
## DomD         0.072  0.141 -0.465       
## Ssnwntr:DmD  0.018 -0.414  0.233 -0.730
print(anova(add, type=2), digits=3)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##            Sum Sq Mean Sq NumDF DenDF F.value  Pr(>F)    
## Season      0.107   0.107     1   112    2.08 0.15217    
## Light       0.422   0.422     1   114    8.22 0.00494 ** 
## Dom         0.308   0.308     1   113    6.01 0.01580 *  
## Season:Dom  0.788   0.788     1   112   15.35 0.00015 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-emmeans(full, ~Season*Dom)
CLD(posthoc, Letters=letters)
##  Season Dom    emmean         SE    df  lower.CL  upper.CL .group
##  winter D   0.2892355 0.08014858  7.79 0.1035443 0.4749266  a    
##  summer D   0.4009953 0.13294032 43.72 0.1330241 0.6689665  ab   
##  summer C   0.4314180 0.06772960  4.06 0.2444074 0.6184287  ab   
##  winter C   0.5608237 0.07554987  6.06 0.3763859 0.7452615   b   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
plot(allEffects(add), ylab="d15N-hs", par.strip.text=list(cex=0.7))

ranef(add)
## $Location
##        (Intercept)
## F1-R46 -0.03524029
## F8-R10 -0.06509519
## HIMB   -0.05945755
## R42     0.15979303
rand(add)
## Analysis of Random effects Table:
##          Chi.sq Chi.DF p.value    
## Location   12.1      1   5e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fixef(add)
##       (Intercept)      Seasonwinter             Light              DomD 
##        0.55502427        0.15503519       -0.01526003        0.09096478 
## Seasonwinter:DomD 
##       -0.38043971
sjp.lmer(add, y.offset = .4)

host C:N

###### host..C.N --
Y<-model.data$host..C.N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     6 -271.49 -254.82 141.75  -283.49                         
## object 10 -266.00 -238.21 143.00  -286.00 2.5084      4     0.6431
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: -252.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8291 -0.7150 -0.1038  0.6546  3.2395 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Location (Intercept) 0.0008875 0.02979 
##  Residual             0.0052739 0.07262 
## Number of obs: 119, groups:  Location, 4
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  1.946e+00  2.324e-02 1.054e+01  83.720 4.44e-16 ***
## Seasonwinter 4.134e-03  1.512e-02 1.139e+02   0.274    0.785    
## Light        1.116e-03  1.649e-03 1.137e+02   0.677    0.500    
## DomD         1.002e-02  1.670e-02 1.141e+02   0.600    0.550    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.539              
## Light       -0.633  0.470       
## DomD         0.135 -0.256 -0.441
print(anova(add, type=2), digits=4)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##           Sum Sq   Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.0003945 0.0003945     1 113.9  0.0748  0.785
## Light  0.0024172 0.0024172     1 113.7  0.4583  0.500
## Dom    0.0018996 0.0018996     1 114.1  0.3602  0.550
plot(allEffects(add), ylab="C:N host", par.strip.text=list(cex=0.7))

ranef(add)
## $Location
##        (Intercept)
## F1-R46 -0.01031414
## F8-R10  0.01985775
## HIMB    0.02386210
## R42    -0.03340572
rand(add)
## Analysis of Random effects Table:
##          Chi.sq Chi.DF p.value   
## Location   6.85      1   0.009 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fixef(add)
##  (Intercept) Seasonwinter        Light         DomD 
##  1.945641081  0.004134462  0.001116126  0.010024231
sjp.lmer(add, y.offset = .4)

symbiont C:N

####### symb..C.N -- 
Y<-model.data$symb..C.N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model
## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
##        Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## ..1     6 -210.60 -193.98 111.30  -222.60                         
## object 10 -205.91 -178.21 112.96  -225.91 3.3076      4     0.5077
summary(add)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
##    Data: model.data
## 
## REML criterion at convergence: -192.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8681 -0.6308 -0.1693  0.6705  2.9795 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Location (Intercept) 4.491e-18 2.119e-09
##  Residual             9.188e-03 9.585e-02
## Number of obs: 118, groups:  Location, 4
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   2.053e+00  2.219e-02  1.140e+02  92.531   <2e-16 ***
## Seasonwinter -2.958e-02  1.956e-02  1.140e+02  -1.513    0.133    
## Light         4.094e-06  1.982e-03  1.140e+02   0.002    0.998    
## DomD          5.402e-03  2.150e-02  1.140e+02   0.251    0.802    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ssnwnt Light 
## Seasonwintr -0.680              
## Light       -0.800  0.426       
## DomD         0.116 -0.226 -0.400
print(anova(add, type=2), digits=5)
## Analysis of Variance Table of type II  with  Satterthwaite 
## approximation for degrees of freedom
##           Sum Sq   Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.0210230 0.0210230     1   114 2.28813 0.1331
## Light  0.0000000 0.0000000     1   114 0.00000 0.9984
## Dom    0.0005801 0.0005801     1   114 0.06314 0.8021
plot(allEffects(full), ylab="C:N symb", par.strip.text=list(cex=0.7))

Figures

Physiology

#########################
# Physiology
df.fig<-data.trim
# figure layout
layout(matrix(c(1,1,2,3), 2,2, byrow=TRUE))
par(mar=c(5,4.5,2,2))

##################
################## 

## season relationship
data.summer<-df.fig[df.fig$Season=="summer", ]
data.winter<-df.fig[df.fig$Season=="winter", ]

## Site relationship
F8.R10.df<-df.fig[df.fig$Location=="F8-R10", ]
F1.R46.df<-df.fig[df.fig$Location=="F1-R46", ]
HIMB.df<-df.fig[df.fig$Location=="HIMB", ]
R42.df<-df.fig[df.fig$Location=="R42", ]

## Symbiont relationship
C.sum.df<-df.fig[(df.fig$Dom=="C" & df.fig$Season=="summer") ,]
C.win.df<-df.fig[(df.fig$Dom=="C" & df.fig$Season=="winter") ,]; C.win.df<-na.omit(C.win.df)
D.sum.df<-df.fig[(df.fig$Dom=="D" & df.fig$Season=="summer") ,]; D.sum.df<-na.omit(D.sum.df)
D.win.df<-df.fig[(df.fig$Dom=="D" & df.fig$Season=="winter") ,]; D.win.df<-na.omit(D.win.df)

##################
# Fig: chlorophyll (total) over season and depth
################## 

plot(chltot~newDepth, data=data.winter,col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")), 
     main="chl: Depth and Season",
     xlim=c(0,11), ylim=c(1,17))
abline((lm(chltot~newDepth, data=data.winter)), col='dodgerblue3', lwd=2)
points(chltot~newDepth, data=data.summer, col="tomato2", pch=16, cex=0.7)
abline((lm(chltot~newDepth, data=data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(data.winter$chltot), ylim=c(0,0.3), xlim=c(0, 18), col="dodgerblue3", main="chl: Seasons", 
     xlab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")))
lines(density(data.summer$chltot), col="tomato2")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.35, -0.1))


###### density plot: by Site
plot(density(F1.R46.df$chltot), ylim=c(0,0.3), xlim=c(0, 18), col="tomato2", main="chl: Sites", 
     xlab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")), ylab="Density")
lines(density(R42.df$chltot), col="skyblue3")
lines(density(F8.R10.df$chltot), col="springgreen4")
lines(density(HIMB.df$chltot), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/physiology/PanKB.chl.pdf", width=9, height=7)



##################
# Fig: chlorophyll (total) over season and depth
################## 

plot(chltot~Light, data=C.sum.df,col="tomato2", pch=16, cex=0.7, 
     xlab="Light (DLI)", ylab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")), 
     main="chl: Light and Season",
     xlim=c(0,20), ylim=c(1,12))
abline((lm(chltot~Light, data=C.sum.df)), col='tomato2', lwd=2)
points(chltot~Light, data=C.win.df, col="dodgerblue3", pch=16, cex=0.7)
abline((lm(chltot~Light, data=C.win.df)), col='dodgerblue3', lwd=2)
points(chltot~Light, data=D.sum.df, col="mediumseagreen", pch=16, cex=0.7)
abline((lm(chltot~Light, data=D.sum.df)), col='mediumseagreen', lwd=2)
points(chltot~Light, data=D.win.df, col="orchid", pch=16, cex=0.7)
abline((lm(chltot~Light, data=D.win.df)), col='orange', lwd=2)
legend("topright", c("C-sum", "C-win", "D-sum", "D-win"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3","mediumseagreen", "orange"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.05))

plot(density(C.sum.df$chltot), col="tomato2", pch=16, cex=0.7, 
     xlab="Light (DLI)", ylab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")), 
     main="chl: Light and Season", ylim=c(0, 0.4), xlim=c(0,15), lwd=2)
lines(density(C.win.df$chltot), col="dodgerblue3", lwd=2)
lines(density(D.sum.df$chltot), col="orange", lwd=2)
lines(density(D.win.df$chltot), col="mediumseagreen", lwd=2)
legend("topright", c("C-sum", "C-win", "D-sum", "D-win"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3", "orange", "mediumseagreen"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.05))


###### density plot: by season
plot(density(data.winter$chltot), ylim=c(0,0.3), xlim=c(0, 18), col="dodgerblue3", main="chl: Seasons", 
     xlab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")))
lines(density(data.summer$chltot), col="tomato2")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.35, -0.1))


###### density plot: by Site
plot(density(F1.R46.df$chltot), ylim=c(0,0.3), xlim=c(0, 18), col="tomato2", main="chl: Sites", 
     xlab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")), ylab="Density")
lines(density(R42.df$chltot), col="skyblue3")
lines(density(F8.R10.df$chltot), col="springgreen4")
lines(density(HIMB.df$chltot), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/physiology/PanKB.chl.pdf", width=9, height=7)





##################
# Fig: chlorophyll/cell over season and depth
##################  
plot(chlcell~newDepth, data=data.winter,col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", 
     ylab=expression(paste("pg chlorophyll cell"^-1, sep="")), main="chl/cell: Depth and Season",
     xlim=c(0,10), ylim=c(0,15))
abline((lm(chlcell~newDepth, data=data.winter)), col='dodgerblue3', lwd=2)
points(chlcell~newDepth, data=data.summer, col="tomato2", pch=16, cex=0.7)
abline((lm(chlcell~newDepth, data=data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(data.winter$chlcell), ylim=c(0,0.4), xlim=c(0, 15), col="dodgerblue3", main="chl/cell: Seasons", xlab=expression(paste("pg chlorophyll cell"^-1, sep="")))
lines(density(data.summer$chlcell), col="tomato2")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(F1.R46.df$chlcell), ylim=c(0,0.4), xlim=c(0,15), col="tomato2", main="chl/cell: Sites", xlab=expression(paste("pg chlorophyll cell"^-1, sep="")))
lines(density(R42.df$chlcell), col="skyblue3")
lines(density(F8.R10.df$chlcell), col="springgreen4")
lines(density(HIMB.df$chlcell), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/physiology/PanKB.chlcell.pdf", width=9, height=7)



##################
# Fig: symbionts over season and depth
##################        

plot((zoox/10^6)~newDepth, data=data.winter,col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab=(expression(paste(~italic("Symbiodinium") ~(10^6~cells~cm^-2), sep=""))), xlim=c(0,10), ylim=c(0, 6),
     main= "zoox: Depth and Season")
abline((lm((zoox/10^6)~newDepth, data=data.winter)), col='dodgerblue3', lwd=2)
points((zoox/10^6)~newDepth, data=data.summer, col="tomato2", pch=16, cex=0.7)
abline((lm((zoox/10^6)~newDepth, data=data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(data.summer$zoox/10^6), col="tomato2", main="zoox: Seasons", 
     xlab=(expression(paste(~italic("Symbiodinium") ~(10^6~cells~cm^-2), sep=""))), xlim=c(0, 7))
lines(density(data.winter$zoox/10^6), col="dodgerblue3")
legend('topright', c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(F1.R46.df$zoox/10^6), col="tomato2", main="zoox: Sites", 
     xlab=(expression(paste(~italic("Symbiodinium") ~(10^6~cells~cm^-2), sep=""))), xlim=c(0, 7))
lines(density(R42.df$zoox/10^6), col="skyblue3")
lines(density(F8.R10.df$zoox/10^6), col="springgreen4")
lines(density(HIMB.df$zoox/10^6), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/physiology/PanKB.zoox.pdf", width=9, height=7)



##################
# Fig: biomass over season and depth
################## 

plot(biomass~newDepth, data=data.winter,col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste("Total biomass", ~(mg~cm^-2), sep="")),
     xlim=c(0,10), ylim=c(5,60), 
     main="biomass: Depth and Season")
abline((lm(biomass~newDepth, data=data.winter)), col='dodgerblue3', lwd=2)
points(biomass~newDepth, data=data.summer, col="tomato2", pch=16, cex=0.7)
abline((lm(biomass~newDepth, data=data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.1))


###### density plot: by season
plot(density(data.winter$biomass), ylim=c(0,0.08), xlim=c(0, 70), col="dodgerblue3", main="biomass: Seasons", xlab=expression(paste("Total biomass", ~(mg~cm^-2), sep="")))
lines(density(data.summer$biomass), col="tomato2")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.35, -0.1))


###### density plot: by Site
plot(density(F1.R46.df$biomass), ylim=c(0,0.1), xlim=c(0,70), col="tomato2", main="biomass: Sites", 
     xlab=expression(paste("Total biomass", ~(mg~cm^-2), sep="")))
lines(density(R42.df$biomass), col="skyblue3")
lines(density(F8.R10.df$biomass), col="springgreen4")
lines(density(HIMB.df$biomass), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), 
       y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/physiology/PanKB.biomass.pdf", width=9, height=7)

Isotopes

### ISOTOPES ##
############### 

##################
# Fig: d13C host over season and depth
##################        

plot(host..d13C~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)")),
     xlim=c(0,10), ylim=c(-20,-10),
     main= "d13C host: Depth and Season")
abline(lm(host..d13C~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(host..d13C~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(host..d13C~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$host..d13C)), col="tomato2", main="d13C-host: Seasons", 
     xlab=expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)")), xlim=c(-20, -8))
lines(density(na.omit(data.winter$host..d13C)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$host..d13C)), col="tomato2", main="d13C-host: Sites", 
     xlab=expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)")), ylim=c(0,0.4), xlim=c(-20, -8))
lines(density(na.omit(R42.df$host..d13C)), col="skyblue3")
lines(density(na.omit(F8.R10.df$host..d13C)), col="springgreen4")
lines(density(na.omit(HIMB.df$host..d13C)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.d13C-host.pdf", width=9, height=7, encod="MacRoman")



##################
# Fig: d13C symb over season and depth
##################        

plot(symb..d13C~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)")),
     xlim=c(0,10), ylim=c(-20,-10),
     main= "d13C symb: Depth and Season")
abline(lm(symb..d13C~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(symb..d13C~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(symb..d13C~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$symb..d13C)), col="tomato2", main="d13C-symb: Seasons", 
     xlab=expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)")), xlim=c(-20, -8))
lines(density(na.omit(data.winter$symb..d13C)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$symb..d13C)), col="tomato2", main="d13C-symb: Sites", 
     xlab=expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)")), ylim=c(0,0.4), xlim=c(-20, -8))
lines(density(na.omit(R42.df$symb..d13C)), col="skyblue3")
lines(density(na.omit(F8.R10.df$symb..d13C)), col="springgreen4")
lines(density(na.omit(HIMB.df$symb..d13C)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.d13C-symb.pdf", width=9, height=7, encod="MacRoman")


##################
# Fig: d13C skeleton over season and depth
##################        

plot(d13C..skel~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste(delta^{13}, C[Skel], " (\u2030, V-PDB)")),
     xlim=c(0,10), ylim=c(-6,2),
     main= "d13C skel: Depth and Season")
abline(lm(d13C..skel~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(d13C..skel~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(d13C..skel~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$d13C..skel)), col="tomato2", main="d13C-skel: Seasons", 
     xlab=expression(paste(delta^{13}, C[Skel], " (\u2030, V-PDB)")), ylim=c(0, 0.5), xlim=c(-8, 6))
lines(density(na.omit(data.winter$d13C..skel)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$d13C..skel)), col="tomato2", main="d13C-skel: Sites", 
     xlab=expression(paste(delta^{13}, C[Skel], " (\u2030, V-PDB)")), ylim=c(0,0.5), xlim=c(-8, 6))
lines(density(na.omit(R42.df$d13C..skel)), col="skyblue3")
lines(density(na.omit(F8.R10.df$d13C..skel)), col="springgreen4")
lines(density(na.omit(HIMB.df$d13C..skel)), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.d13C-skel.pdf", width=9, height=7, encod="MacRoman")


##################
# Fig: d13C host-symb over season and depth
##################        

plot(d13C..host.sym~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)")), 
     xlim=c(0,10), ylim=c(-2.5,3),
     main= "d13C h-s: Depth and Season")
abline(lm(d13C..host.sym~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(d13C..host.sym~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(d13C..host.sym~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$d13C..host.sym)), col="tomato2", main="d13C h-s: Seasons", 
     xlab=expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)")), ylim=c(0,1.5), xlim=c(-3, 5))
lines(density(na.omit(data.winter$d13C..host.sym)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$d13C..host.sym)), col="tomato2", main="d13C h-s: Sites", 
     expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)")), ylim=c(0,1.7), xlim=c(-3,5))
lines(density(na.omit(R42.df$d13C..host.sym)), col="skyblue3")
lines(density(na.omit(F8.R10.df$d13C..host.sym)), col="springgreen4")
lines(density(na.omit(HIMB.df$d13C..host.sym)), col="purple")
legend("topright",c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.d13Ch-s.pdf", width=9, height=7, encod="MacRoman")




##################
# Fig: d15N host over season and depth
##################        

plot(host..d15N~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste(delta^{15}, N[H], " (\u2030, air)")), xlim=c(0,10), ylim=c(2,8),
     main= "d15N host: Depth and Season")
abline(lm(host..d15N~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(host..d15N~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(host..d15N~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$host..d15N)), col="tomato2", main="d15N-host: Seasons", 
     xlab=expression(paste(delta^{15}, N[H], " (\u2030, air)")), xlim=c(0, 12), ylim=c(0, 1))
lines(density(na.omit(data.winter$host..d15N)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$host..d15N)), col="tomato2", main="d15N-host: Sites", xlab=expression(paste(delta^{15}, N[H], " (\u2030, air)")), ylim=c(0,2), xlim=c(0, 12))
lines(density(na.omit(R42.df$host..d15N)), col="skyblue3")
lines(density(na.omit(F8.R10.df$host..d15N)), col="springgreen4")
lines(density(na.omit(HIMB.df$host..d15N)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.d15N-host.pdf", width=9, height=7, encod="MacRoman")


##################
# Fig: d15N symb over season and depth
##################        

plot(symb..d15N~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab =expression(paste(delta^{15}, N[S], " (\u2030, air)")),
     xlim=c(0,10), ylim=c(2,7),
     main= "d15N symb: Depth and Season")
abline(lm(symb..d15N~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(symb..d15N~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(symb..d15N~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$symb..d15N)), col="tomato2", main="d15N-symb: Seasons", 
     xlab= expression(paste(delta^{15}, N[S], " (\u2030, air)")), xlim=c(0, 9), ylim=c(0, 0.9))
lines(density(na.omit(data.winter$symb..d15N)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$symb..d15N)), col="tomato2", main="d15N-symb: Sites", 
     xlab= expression(paste(delta^{15}, N[S], " (\u2030, air)")),  xlim=c(0, 9), ylim=c(0,1.8))
lines(density(na.omit(R42.df$symb..d15N)), col="skyblue3")
lines(density(na.omit(F8.R10.df$symb..d15N)), col="springgreen4")
lines(density(na.omit(HIMB.df$symb..d15N)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.d15N-symb.pdf", width=9, height=7, encod="MacRoman")



##################
# Fig: d15N host.symb over season and depth
##################        

plot(d15N..host.sym~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste(delta^{15}, N[H-S], " (\u2030, air)")),
     xlim=c(0,10), ylim=c(-1,2),
     main= "d15N h-s: Depth and Season")
abline(lm(d15N..host.sym~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(d15N..host.sym~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(d15N..host.sym~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$d15N..host.sym)), col="tomato2", main="d15N h-s: Seasons", 
     xlab=expression(paste(delta^{15}, N[H-S], " (\u2030, air)")), xlim=c(-1, 3), ylim=c(0,2))
lines(density(na.omit(data.winter$d15N..host.sym)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.3, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$d15N..host.sym)), col="tomato2", main="d15N h-s: Sites", xlab=expression(paste(delta^{15}, N[H-S], " (\u2030, air)")), ylim=c(0,2), xlim=c(-1,3))
lines(density(na.omit(R42.df$d15N..host.sym)), col="skyblue3")
lines(density(na.omit(F8.R10.df$d15N..host.sym)), col="springgreen4")
lines(density(na.omit(HIMB.df$d15N..host.sym)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.d15Nh-s.pdf", width=9, height=7, encod="MacRoman")



##################
# Fig: C.N host season and depth
##################       
plot(host..C.N~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste("C:N"[H])),
     xlim=c(0,10), ylim=c(5,10),
     main= "C:N-host Depth and Season")
abline(lm(host..C.N~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(host..C.N~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(host..C.N~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$host..C.N)), col="tomato2", main="C:N-host Depth and Season", 
     xlab=expression(paste("C:N"[H])), ylim=c(0, 1), xlim=c(4, 10))
lines(density(na.omit(data.winter$host..C.N)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$host..C.N)), col="tomato2", main="C:N-host Depth and Season", 
     xlab=expression(paste("C:N"[H])), ylim=c(0,1.5), xlim=c(4, 10))
lines(density(na.omit(R42.df$host..C.N)), col="skyblue3")
lines(density(na.omit(F8.R10.df$host..C.N)), col="springgreen4")
lines(density(na.omit(HIMB.df$host..C.N)), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.C.Nhost.pdf", width=9, height=7, encod="MacRoman")



##################
# Fig: C.N symb season and depth
##################       
plot(symb..C.N~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7, 
     xlab="Depth (m)", ylab = expression(paste("C:N"[S])),
     xlim=c(0,10), ylim=c(5,14),
     main= "C:N-symb Depth and Season")
abline(lm(symb..C.N~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(symb..C.N~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(symb..C.N~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))

###### density plot: by season
plot(density(na.omit(data.summer$symb..C.N)), col="tomato2", main="C:N-symb Depth and Season", 
     xlab=expression(paste("C:N"[S])), ylim=c(0, 0.8), xlim=c(5, 15))
lines(density(na.omit(data.winter$symb..C.N)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2), 
       col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))

###### density plot: by Site
plot(density(na.omit(F1.R46.df$symb..C.N)), col="tomato2", main="C:N-symb Depth and Season", 
     xlab=expression(paste("C:N"[S])), ylim=c(0,0.8), xlim=c(5, 15))
lines(density(na.omit(R42.df$symb..C.N)), col="skyblue3")
lines(density(na.omit(F8.R10.df$symb..C.N)), col="springgreen4")
lines(density(na.omit(HIMB.df$symb..C.N)), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
       col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))

dev.copy(pdf, "figures/isotope/PanKB.C.Nsym.pdf", width=9, height=7, encod="MacRoman")